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A MESH GRADIENT TECHNIQUE FOR 


NUMERICAL OPTIMIZATION 


ABSTRACT 

by 

EDWARD ALLEN WILLIS, JR. 


This paper deals with a class of successive-improvement optimi- 
zation methods in which directions of descent are defined in the 
state space along each trial trajectory. The given problem is first 
decomposed into two discrete levels by imposing mesh points. Level I 
then consists of running optimal subarcs between each successive 
pair of mesh points. For normal systems, these optimal two-point 
boundary value problems can be solved by following a routine pre- 
scription if the mesh spacing is sufficiently close . A spacing cri- 
terion is given. Under appropriate conditions, the criterion value 
depends only on the coordinates of the mesh points, and its gradient 
with respect to those coordinates may be defined by interpreting the 
adjoint variables as partial derivatives of the criterion value 
function. In Level II , the gradient data is used to generate im- 
provement steps or search directions in the state space which satisfy 
the boundary values and constraints of the given problem. The family 
of feasible varied trajectories thus constructed converges to the 
"nearest" locally-optimum trajectory, if any such exist. 
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1. SUMMARY 

This paper deals with a class of successive-improvement optimiza- 
tion methods in which directions of descent are defined iii the state 
space along each trial trajectory. The given problem is first decom- 
posed into two discrete levels by imposing mesh points, Level I then 
consists of running optimal subarcs between each successive pair of 
mesh points. For normal systems, these optimal two-point boundary 
value problems can be solved by following a routine prescription if 
the mesh spacing is sufficiently close , A spacing criterion is 
given. Under appropriate conditions, the criterion value depends 
only on the coordinates of the mesh points, and its gradient with 
respect to those coordinates may be defined by interpreting the ad- 
joint variables as partial derivatives of the criterion value func- 
tion. In Level II , the gradient data is used to generate improvement 
steps or search directions in the state space which satisfy the 
boundary values and constraints of the given problem. The family of 
feasible varied- trajectories thus constructed converges to the "near- 
est" locally-optimum trajectory, if any such exist. 

This approach reduces the typical deterministic optimal control 
problem to an exercise in the classical theory of maxima and minima. 

It also leads to a class of apparently novel algorithms for optimal 
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trajectory computations. These differ as to their ultimate rates of 
convergence and computational sophistication, ljut all share the same 
basic operation “ namely, decomposition and variation of the state 
trajectory itself* This derives maximum benefit from the initially- 
given constraint and boundary value information* It is also useful 
for solving computationally-unstable problems (where the integration 
interval is long compared to state or adjoint system time constants) , 
by taking the mpsh-spacing small enough* Therefore, the present 
class of methods prove relatively effective when applied tp unstable 
problems or problems wifh numerous boundary vplues and constraints , 
as illustrated by several examples. 



2. INTRODUCTION 


Many important physical processes are describable by systems of 
deterministic ordinary differential equations, Since there are often 
more variables than equations, it makes sense to use some of the ex- 
tra degrees of freedom as control variables - in order to optimize a 
criterion of merit and to meet the boundary values prescribed for the 
process. Thus (with many technicalities deferred until the formal 
"statement of the problem" in Section 2,1) the systems studied may 
typically have the form 


minimize J = 



f Q (x(t) ,u(t) , t)dt 


with 

x(t) = f(x(t) ,u(t) ,t) 

subject to various boundary conditions and constraints on the state 
vector x(t) and the control vector u(t). 

From a mathematical viewpoint, modern theories, spch as, L. S. 
Pontryagin's Maximum Principle (ref. 137; references are listed 
alphabetically in Chapter 6, Bibliography) give enough "necessary 
conditions" to define solutions to most practical problems. Unfor- 


tunately, it is not always easy to implement these conditions in 
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practical numerical studies. Experience indicates that the difficul- 
ties experienced in conducting numerical studies are usually due to 
one or both of the following causes: 

(a) Numerical instability results, when the integration interval 
exceeds the dominant state or adjoint time constant, in computer over- 
flows, inordinate sensitivity of final condition to small perturba- 
tion of the initial data, etc. That this is a widespread problem may 
be inferred by recalling that for a simple linear system, each stable 
state equation will give rise to an unstable adjoint equation, and 
conversely. 

(b) Numerous constraints on the state trajectory (c.f . fig. 

2-1), involving point or path, equality or inequality relations, nor- 
mally entail the enforcement of additional relations, such as the 
"jump" and "corner" conditions or appropriate forms of the transver- 
sality conditions. Typically, one must introduce and determine an 
extra set of auxiliary multipliers for every intermediate constraint, 
thus adding to the computational burden. Philosophically, one would 
expect the presence of intermediate constraints to help, not 

hinder the optimization process by reducing the range and dimension 
of the space to be searched. 

These factors, which are greatly compounded by the nonlinearity 
and high dimensionality typical of practical problems, are respon- 
sible for phenomena such as numerical instability and multiple or 
nonexistent solutions which seriously hinder the effective conduct 








and p € is a vector of constant (but adjustable) design 
parameters . 

In addition there are terminal- and intermediate-point equality 
constraints, i.e., 

^V-^^’t) (4) 

where tg < t^ ... t R = t f and the are dimensional 

smooth mani folds (0 _< n^ <. N) . 



7 


Finally, state-variable inequality constraints may also apply, 


i.e. 


x(t) £ int T 


(5) 


where F represents one or more excluded regions of state 

space. 

2.2 Existing Methods of Solution 

Necessary conditions and sufficient .conditions for basic forms 
of the problem are well in hand, as may be judged from the appear- 
ance recently of comprehensive texts such as Athans & Falb (ref. 4), 
Bryson & Ho (ref. 24), and Lee & Markus (ref. 108). Several basic 
computational approaches are presently available for producing num- 
erical solutions in specific cases. Since recent surveys (e.g., 
refs. 3, 54, 96, 132 and 149) display literally hundreds of individ- 
ual contributions, the present discussion will be limited to broad 
categories, with details left to the indicated references.* 

Conceptually, the most powerful method now available is Dynamic 
Programming. Although a valuable theoretical tool, its computational 
application have been very limited due to extremely high storage re- 
quirements. Hence, despite the appearance of D.P. algorithms with de- 
creased storage requirements (refs. 97 and 98), it is felt that iter- 
ative successive-improvement are more promising as a class. 

First-order schemes such as Kelley’s control iteration (ref. 81) 
are perhaps useful for generating initial guesses for higher ordered 
schemes. Unfortunately, they converge slowly in the neighborhood of 
the solution and tend to produce inaccurate adjoint trajectories. 

* 

References are listed and categorized in Chapter 6, Bibliography. 
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These disadvantages, which are pointed out and discussed in connec- 
tion with Figures 5 and 6 of reference 100, render first-order tech- 
niques generally unsuitable for the present needs, and they will not 
be further discussed. 

Newton-type methods achieve quadratic convergence in the neigh- 
borhood of the true solution, by making use of linear system theory. 
The Transition Matrix approach (c.f. ref. 24) is one example of a 
Newtonian (second order) boundary iteration technique. The major 
disadvantage of this method is its tendency toward numerical instab- 
ility. As previously explained either the state or adjoint equations 
are necessarily unstable. This can result in the state and adjoint 
variables being of entirely different orders of magnitude (e.g., e at 

“clt 

vs e ) - which implies the transition matrix is ill-conditioned - 
when the time interval significantly exceeds the system's dominant 
time constant. 

To reduce numerical-instability problems, the "backward-sweep" 
approach can be applied to both boundary-iteration (ref. 24) and 
control-iteration (ref. 123) techniques. As applied to boundary 
iteration, this entails integrating the state, adjoint and Ricatti 
matrix equations backwards from assumed terminal conditions. Term- 
inal values are then adjusted by iteration so that the desired ini- 
tial conditions are achieved. The advantage of this approach is 
that the adjoint equations are being integrated in the "stable" di- 
rection. The Ricatti equation may also be tested for system "nor- 
mality" and the presence of conjugate points (c.f. Section 8.2, 
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Appendix). On the other hand, even though the stability problem may 
be somewhat eased in some cases, it is by no means eliminated. This 
is because the canonical equations necessarily contain both stable 
and unstable components and these are merely interchanged by revers- 
ing the sense of integration. 

When control rather than boundary-value iteration is used, the 
state equations are integrated forward using an assumed control his- 
tory and the resulting trajectory stored; the adjoint and Ricatti 
matrix equation are then integrated backwards to determine an im- 
proved control history which is both more nearly optimal and more 
nearly feasible. Stability is greatly improved because the state 
and adjoint equations are decoupled during any one pass. Neverthe- 
less, a fairly good initial guess is still required for the entire 
control history, and (short of first making use of a lower-order 
scheme) there is no obvious or systematic waj T to obtain one. In ad- 
dition, this technique requires both forward and backward integration 
passes with storage of the entire state and adjoint trajectories and 
Ricatti matrix histories. It requires fixed-step-size integration 
schemes which, in general, are less efficient than variable step size 
schemes. Finally, the process does not inherently meet all boundary 
conditions and an additional multiplier must be determined for each 
"hard" terminal constraint (i • e • > one which must be satisfied ex- 
actly). On balance, this approach is rather attractive for free- 
terminal problems but rapidly loses it appeal when "hard" con- 


straints are added. 
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Finally, the best known of the Newtonian functional approxima- 
tion techniques is that of Quasilinearization (ref. 14). The proce- 
dure is to simply linearize the canonical equations around trial 
state/adjoint trajectories. The following linear perturbation equa- 


tions 


<5x l(x.4o - 1 

l r i l 

<5^ [g(x j .,^ i ) - ip ± 



* 1+1 (t) = *i (t) + ^i+i (t) = ^i (t) + (8) 

are solved itteratively until a suitable norm, e. g., max' 6x (t) , 

t€^ 

becomes "sufficiently" small. The procedure unfortunately is subject 
to instability whenever the integration interval significantly ex- 
ceeds the dominant time constant and also is dependent on having a 

good initial guess for both the state and adjoint trajectories. The 

/ 

technique also requires storage of the entire state-adjoint trajec- 
tories and is not adaptable to efficient, variable-step-size numeri- 
cal integration routines. 

From the preceding discussion it may be concluded that none of 
the existing methods seems to be specifically addressed to the most 
general, unstable, and highly constrained ,form of the problem. More- 
over, when nonlinear effects are added to the previously-mentioned 
instability and constraint problems (c.f. Fig. 2-1), it is readily 
understandable why computational efficiency is still a major consid- 
eration despite the capabilities of modern computer facilities. 
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State variable inequality constraints or intermediate-equality con- 
straints are particularly troublesome in that they entail additional 
mathematical conditions and may require a laborious "patching to- 
gether" of free and constrained arcs. Philosophically, it would 
seem that the additional a priori information represented by the 
state constraints should be usable to help , not hinder, in computing 
the solution. 



Based on the preceding remarks , the present approach was de- 
signed with two primary objectives in view: (a) to take maximum ad- 

vantage of the a priori data given about the problem in the form of 
intermediate state-variable constraints; and (b) to eliminate or 
avoid the problem of numerical instability. 

In essence, the present approach relies upon the principle 
of decomposition to divide an intractable given problem into nested 


levels of individually-manageable subproblems. Specifically the pro- 
cedure is to impose a set of mesh points {x^;t^} leading from (Xq^q) 
to (x K ;t K ), certain members of which belong to the constraint mani- 
folds^. This resolves the problem into two levels- 


Level I : Select an initial mesh, { (x 0 ;t Q ) , (x-^t^) •••• (x^t^l 

and (with the mesh-points held fixed) solve the resulting sequence 
of Two-Point Boundary Value Problems (TPBVP) using any suitable nu- 
merical optimization technique. These solutions may be readily 
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computed by taking closely-spaced mesh points; this avoids the sta- 
bility problems often associated with initial-value methods, and re- 
sults in an initial feasible trajectory satisfying the given phys- 
ical boundary values and constraints. 

Level II : Optimize the mesh-point locations so as to minimize 

J , using any suitable mathematical programming technique. Note that 
here the physical constraints and boundary values merely reduce the 
range and/or dimensionality of the search. 

2.4 Comments 

(a) If a certain controllability condition is satisfied, numeri- 
cal instabilities may always be avoided in Level I by taking a suffi- 
ciently close mesh spacing. This is an important point, because 
Level I is nested inside Level II and must operate through many 
cycles for every step in the Level II search. Thus, it is essential 
for the present method to have a reliable and efficient method of 
solving short subarc TPBVP's. 

(b) If the TPBVP solutions are optimal, then J is a function 
of the mesh-point coordinates only. The gradient of J may be read- 
ily defined in terms of the adjoint-variable and Hamiltonian-f unction 
discontinuities at each mesh point. Furthermore, the Hessian 
matrix may be derived from the subarc transition matrices. Thus, 
Level II consists of a conventional, well-posed mathematical program- 
ming problem for which a very satisfactory and complete theory is 
available - c.f. Fiacco & McCormick (ref. 45). 



(c) The selection of an initial set of mesh points is to some 
extent a matter of judgment. As a minimum, the initial and final 
points, plus one point for every junction with a constraint mani- 
fold, must be included. Additional unconstrained points may be in- 
serted to satisfy the mesh-spacing criterion for stability (to be 
developed later) . 

(d) The tractability of closely spaced TPBVP's depends on a 
controllability assumption - namely, that each mesh point is in fact 
reachable from the preceding one. If this is not true in general 
(as is apt to be the case with fixed-time, bounded-control prob- 
lems) , then care must be used to avoid the appearance of abnormal 
(i.e., unfeasible) subarcs. In many cases, the TPBVP's may be re- 
formulated in terms of quantities that are attainable with the avail- 
able control. This would typically involve dropping redundant or un- 
controllable coordinates, substituting iteration parameters, etc. 

See reference 115 for instance. Otherwise, penalty functions may be 
used to weigh temporarily-unavoidable boundary value or control- 
constraint violations. 

(e) The Level I and Level II operations individually make use 

of established techniques and theories. The novelty and contribution 
of the Mesh-Gradient approach consists in combining these techniques 
in such a way as to avoid numerical instabilities and to make effec- 
tive use of the constraint data specified for the problem. Whereas, 
prior techniques may be broadly classified as involving boundary 



iteration, control iteration or successive approximation (i.e., quas- 
ilinearization) , the characteristic feature of the present approach 
is iteration in the state space . 

(f) This dissertation is organized in the following manner. 

The basic Mesh-Gradient Technique is presented in Chapter 3 (with 
all but the most essential of the supporting material deferred to 
Chapters 8 and 9). Four computational examples are treated in de- 
tail in Chapter 4, while conclusions and recommendations for further 
study are given in Chapter 5. References are listed arid Categorized 
in Chapter 6, main symbols are presented in Chapter 7, and the re-^ 
maining Chapters are Appendices. 



3. THE MESH-GRADIENT TECHNIQUE 

The major steps comprising the "Mesh-Gradient" method are pre- 
sented in this chapter. For the sake of clarity, only its character- 
istic features and underlying assumptions are discussed here, with 
many derivations and other lengthy details deferred until Chapters 8 
and 9. 

As applied to Problem 2.1, the' technique begins by introducing 

a set of mesh points which in some sense lead from the 

given initial conditions to the desired final conditions. The 

points are ordered, so that t A < t, < • • • • t, „ < t, «••• < t = t 

U 1 k-1 k K f 

and initially arranged to lie along the analyst’s best guess at the 
optimal trajectory. The collection of admissible mesh points in 
is denoted by ; other main symbols are listed in Chapter 7. One 
member of belongs to each constraint manifold, and additional, 

unconstrained points may be introduced for computational convenience. 
This defines a sequence of two point boundary value problems (TPBVP) 
connecting adjacent pairs of mesh points. The objective of Level I 
is that each TPBVP solution will be optimal . When this holds true, 
it follows that the criterion value J is a function of the mesh 
point coordinates and times only, and its gradient with respect to 
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these coordinates is well-defined. The notation V^( ) will denote 
the gradient with respect to feasible direction in 3C . Thus, the 
objective of Level II is to achieve that J is minimal with re- 
spect to feasible variations in X , i.e., that V_J = 0 and that 
r pp X 

W 1 is positive definite. 

3 . 1 Level I - Subarc Solutions 

As mentioned in Chapter 2 it is crucial for the present tech- 
nique to have an efficient, reliable, hopefully almost fool-proof 
method of solving TPBVP's. This is because Level I is nested inside 
Level II and must operate through a full cycle (i.e., solve K 
TPBVP's) for every "function-evaluation" in the Level II search. In 
this section we will discuss the chosen technique, and subsequently 
the major conditions and restrictions that apply to it. 

For smooth, well-behaved functions of practical engineering in- 
terest (c.f., the standard definitions and assumptions listed in 
Section 8.1.1) it is well known that an optimal control and its re- 
sponse must satisfy the Maximum Principle (c.f., Sections 8.1.2 and 
8.1.3), an appropriate form of the transversality condition (Section 
3.2.1) and also the convexity or strengthened Legendre-Clebsch con- 
dition, the normality condition and the Jacobi or No-Conjugate- 
Point condition (c.f., Sections 8 .2 . 1-8.2 .3) . From henceforth it 
is assumed that these conditions apply to individual subarcs. 

3.1.1 The Transition Matrix Algorithm 

While many numerical techniques are available for solving two 
point boundary value problems, Newtonian iteration based upon linear 
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perturbation theory (the transition matrix algorithm or TMA) is 
used here because of its efficiency and because the associated 


second-partial matrices can also be used in the upper level calcu- 


lations. Specifically, after introducing the usual adjoint vari- 
ables (4>(t)) and Hamiltonian and ,-using the maximum principle to 

* 

eliminate the control variables, the subarc problem may be expressed 


in the canonical form 


where the 


x(t) = - 20 - = 

'P(t) = - ^0 = g(x(t) J(t) ,g,t) 

dX 

; <V e *? k J «W 4 i k+1 € % +1 


•K. 




a) 


x. €. X , with k = 1,2,"-*"] 

K PP 

With assumed initial values of i|j(t^), Eqs. (1) may be integra- 
ted forward to obtain end points x^(t^ +1 ) (see Fig. 3-1) which in 
general differ from the desired points x^ + ^ , and contributions J k 
to the criterion value. Relative to the given initial points x^ 
and the actually obtained final points x^(t^ + ^) , each J k is at 


least stationary in F 


Moreover, it will be shown that J v de- 


sa -- J k 

pends upon its end-point coordinates only, and that (with suitable 


normalization) the partial derivatives of the integral contribution 

to J k with respect to x(t^) and x *(t k+1 ) are given by 
* 

If is unbounded, it is generally possible (under the hypotheses 
of Section 8.1.1) to solve for u as an explicit, differentiable 
function of x, y, £ and t. If Q is bounded, however, u may 
depend implicitly on some components of ip , i.e., via the location 
of switching boundaries. In that case, some elements of the Jacob- 
ian matrices in Eqs. (6) and (7) below will be represented by im- 
pulse functions , and special techniques must be used to integrate 
the transition matrix Eqs. (4). 



State variable, 


- is - 




FIGURE 3-1. - INITIAL AND FINAL STAGES IN THE LOWER 
LEVEL CYCLE. 
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3x(t fc ) * 


^* (t k+l) 


*0 


where superscripts + and - denote left and right hand limits, re- 


spectively. When P.V.C. terms are present in J, Eq. (2) must be ex- 


tended to include 


respectively. 


-r— (PVC) k and “ — ~ — 
3x ( t k ) 3x * ( Vi) 


(PVC), 


The required values of which cause x^(tj.^) to coincide 

with the desired point x^^, are found by linearizing Eqs. (1) 

around the trajectory ensuing from a previous (j ) estimate ifi„(t,). 

J k 

In the customary way this yields the transition matrix and perturba- 


tion system. 


+ 

fc k+i ' 



! fc k+l 


Arc k contributes 


an amount J, to the criterion value 
k 

For the present purposes it will be convenient to partition the 
transition matrix $(t n ;t) into A N x N blocks, i.e.. 


®(t Q ;t) = 


(3) 
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With 


_d 

dt 


A 

B 


01 

6 


A 

B 

C 

D 

t A ;t 

Y 

6 

t 

C 

D 


( 4 ) 


and 


A B 


C DJ 

- — I f* * f 


= I (the 2N x 2N identity matrix) 


(5) 


'O’ 0 


where a , 3, y, and 6 are NxN Jacobian matrices formed by differ- 

entiating f and g with respect to x and \p in turn, namely 

3^ 0 3? 3g x 3g . 

3x 3 t^ 9x dip 

Recalling the canonical definitions of 1 and g (c.f., Section 
8.1.2), the Jacobians may be expressed (component-wise) as 


= 

a ii 3x. 3 t1> . 

i J 


e 


= 

ij “ dlp^dlp^ 


> 


d 2 *& . = _ 3 2 ^ 

^ij dx_.dx. ij 9^.9x. 


i 3 


i J 


(7) 


Because of the continuity and differentiability of the original func- 
tion ?(x,S,S,t), it may be observed that 3 and y are symmetric, 
i.e. , 


6=6 


Y = Y 


and also that 




6 = - a 
.2 




( 8 ) 


Thus, the Jacobians entail 2N + N independent quantities rather 
than 4N 2 . 
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Then, computing $(tQ;t) as usual and adopting the notation 
that A(t^;t^ + ^) = A^, etc., the solution of the perturbation system 
may be written as 


' : * ( W 


iteration 

j 

,-l 


\ B k 
C k \ 


<Sx(t k ) 

«S^(t k ) 


(9) 


Finally, assuming that ^ exists, we may solve for 6i|i(t k ) with 


<5x(t k ) = 0, yielding an increment 


’ B k lii j ( W 


( 10 ) 


which (for linear dynamics) would zero the predicted terminal error 

6x j+l (t k+l ) - 
3.1.2 Comments 

(1) It is shown in Section 8.3 that the sequence 
(tfc) arbitrary 


Vl (t k> ‘ V*’ + 


(ID 


converges either quadratically or not at all to the solution of (1) . 

(2) Moreover, the assumptions of Section 8.1.1 imply that except 
in the presence of an abnormal or conjugate point trajectory, conver- 
gence can always be obtained if adjacent partition points are prop- 
erly arranged and sufficiently close together (briefly, this means 
that each point is to be attainable from the preceding one and that 
the interval [t^t^^] does not much exceed the dominant local time- 
constant of the canonical system); 
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(3) The results of section 3.2.1 below imply that, under the 


normalization implied in Eq. (2), the blocks of the transition matrix 


$ = 


A B 
C D 


may be interpreted as 


A k (t) ■ 


8i!(t k+l ) 

9x(t k ) 

8x ( t ) 
B k (t) - V - 

9^(t k ) 


C k = 


C k = 


9x(t k ) 

3 » ( W 

4(t k > 


(12) 


for all t € 2T. 

(4) The existence and proper interpretation of the partial de- 
rivatives indicated in this section and below depend upon (a) the 
continuity, differentiability and convexity properties listed in 
Section 8.1.1, and (b) the non-singularity of B k » The latter does 
not necessarily follow from the former, however, and the conditions 
under which B^ or its equivalent exists are further discussed in 
the next section. 

To summarize, the Transition Matrix Algorithm or TMA comprises 
the following steps: 

(a) ^ 0 (t k ) arbitrary 

(b) solve Eqs. (1) to (9) 
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(c) Compute the linear increment Sip . ( t.) and correct ^„(t.) 

J h 3 k 

per Eqs. (10) and (11). 

(d) Iterate steps (a)-(c) until a suitable error norm, say 

5S i (t k+i ) Q6l j (t k+i> (12) 

(where Q is a positive definite matrix) becomes suf- 

I 

ficiently small. 

For a TPBVP, this scheme requires the forward integration of the 

2N cononical equation plus the 4 sets of N x N matrix equations to 

define A, B, C and D. It may be recalled from Eqs. (7) that 

2 

is symmetric so that only 2N + N terms must be computed. 

There is no need to store the trajectory histories. At the final 

time the N x N matrix B must be inverted, and used to compute the 

derived matrices E, F, G and H which will be defined in Section 

3.2. As shown in the Appendix (8 . 2 .4)these matrices collectively con- 
2 

tain only 2N + N independent components. E, F, G, and H may be 
stored in the same locations used for A, B, C, and D, and used 
prior to the next iteration to compute quantities that are of inter- 
est for Level II. 

3.1.3 Conditions For Feasible Solutions 

As has been previously implied, even the generous hypotheses of 
Appendix A are not sufficient to guarantee the existence of a matrix 
B^ which solves Eq. (1) in 3.1.1. This is because, even though 
Problem 2.1 is assumed to be well-posed, i.e., its solution exists 
in principle, the system is in general only partically controllable 
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in the sense of Section 8.1.1. In such a case we have no way of 
knowing apriori whether one mesh point lies in the controllable 
subspace of the next, etc. Thus, the concept of and conditions for 
feasibility are of central importance here. 

A feasible trajectory is one satisfying all specifications of the 
problem except that the criterion J is not necessarily minimized. 
That is, using the attainable set notation defined in Section 8.1.1, a 
trajectory x(t), joining the points x(t k > = ^ and x(t k+1 > = 
x k+ ^, does not exist unless 

^k+i c ^VVW 

or equivalently, 

\ e ^ (t k ;t k+l’\+i> (1) 


An analogous statement with 
tween every pair of distinct 
tory, i.e., the relation 

x(t) ^ 

holds, for any t e. ( t fc , t fc+ ^ 


respect to sets of traversal holds be- 
points belonging to a feasible trajec- 

t k ;t) n ^ ^ t;t k+i»\+i^ ^ 

). 


Local feasibility hypothesis . - The feasibility condition of 
Eq. (2) is obviously satisfied by any pair of points belonging to an 
optimal trajectory. The condition also applies, by continuity, in 
some neighborhood of the optimal trajectory. The local feasibility 
hypothesis consists in assuming that condition (2) applies along an 
arbitrary trial trajectory. This implies that all of the necessary 
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2-point transfers exist. If the hypothesis is untrue for the orig- 
inal problem (as evidenced by the appearance of singular B) , it can 
be relaxed into a controllable version of the problem by temporarily 
introducing point-discontinuity penalty functions. That is, the 
"hard" boundary conditions such as = ^ are removed, and 

in their places a P.V.C. such as the quadratic form 

x “>■ .. T * r - ^ / . ( 3) 

p j “ Vl 1 5 Ix(t k+1) - Vl 1 

(rtiere Q is a positive definite N x N matrix and p^ > 0 is a 
scalar penalty factor) is added to the previous criterion value. 

Then as shown in reference 45, the sequence p^ arbitrary > 0, 
p. > p._^, with lim p. = ", will lead to a trajectory with vanish- 
ingly small mesh point discontinuities, from which the solution to 
the original problem may be recovered. 

In many cases, however, it is possible to avoid the use of pen- 
alty functions by approximately choosing the iteration parameters. 
Also, the adjustable mesh-points should be specified in a manner 
compatible with the initial boundary conditions. For example, for 
time-open or time-optimal problems, the mesh points should not be 
fixed in time. That is, one must formulate the TPBVP's in a fashion 
that makes sense in terms of the given system and the capabilities 
of the control system, by describing the mesh points in terms of 
quantities that are attainable for this particular system . Unfor- 
tunately, this process does require some judgement on the part of 
the analyst and cannot be reduced to a routine prescription. The 
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following considerations, however, are helpful in identifying the 
proper formulation. 

Generalized iteration parameters . - As has been mentioned the 
two point boundary value problems are solved numerically by itera- 
tion. In general, the N terminal conditions are to be met by 
proper choice of N initial conditions. Let the vector d^ ab- 
stractly symbolize the desired terminal condition and p^_^ t ^ ie ^ n_ 
itial condition to be determined. 

Typically, the components of d(t^) are determined by n^ 
terminal equality constraints and N- n^ transversality conditions. 
For example, if requires the first n^ components of x(t^) 

to have fixed final values, the other N- n^ being free, then d(t^) 
has the form 

i(t k > - 


— 


X l <t k ) ’ X 2 (t k ) 


••v (t t >; V+iV 


'W 


(4) 


and must satisfy the terminal condition 


I(t k ) 


*kl ,"*kl ,""*k 

; 0* • • -0 

1 2 

\ 


(5) 


Similarily, represents the initial value of those components 

of the state and adjoint vector which are not defined by boundary 

values or transversality conditions at time t, That is , p, 

k-1 r k-l 

represents N independent initial conditions which lead, via inte- 
gration of the canonical equations, to N dependent final values 
dCt^). Thus, the iteration problem is to find the N-vector 
which will transfer the canonical state-adjoint system to a terminal 
state that satisfies Eqs. (4) and (5). For such a transfer to 
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exist, it is necessary that the symbolic, N x N Jacobian matrix 

S 3(t k ) 


( 6 ) 


8p 


k-1 


be nonsingular. This requirement dictates the choice of the itera- 
tion parameters p^_^. 

It is desirable, whenever possible, to use the initial adjoint 
variables as the p^_^'s; however, other parameters can be used 
when appropriate to avoid abnormal or unfeasible subarcs. For ex- 
ample, the elapsed or final time may be included if the control ef- 
fort available is bounded. This is because an arbitrary mesh point 
x^, even if it satisfies the local feasibility hypothesis of Eq. (2) 
with respect to the preceding and subsequent mesh points, then may 
not belong to «. er (x^_ 1 , t^ t^) for all values of t^. 

For example, even the elementary problem 


with 


minimize J 



u(t) 


,2 

dt 


x(t) = u x(0) = 0 x (1 ) = 1 + €1 
(where 61 > 0) exhibits this behavior if the control magnitude is 
constrained, such that |u(t)| <. 1. No solution then exists unless 

the final time is relaxed at least to 1 + £. 

In more general examples, it may be possible to use components 
of the design parameter vector) "g or of the state or control vec- 
tors, in addition to the final time, to replace some components of 
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ijj (t^). See reference 115 for a relatively thorough discussion of 
this possibility. 

Finally, note that it is possible to investigate the existence 
of B \ at least over a short interval, by examining the transition- 
matrix differential Eq. (4) in 3.1.1. I.e., to the first order in 
time, 

B(t - t Q ) » [a(t Q )B(t 0 ) + 0(t o )D(t o )](t - t Q ) = e(t Q )(t - t Q ) 

(7) 

Hence, if the Jacobian matrix 3(tg) itself is singular, one need 
not solve Eqs. 3. 1.1(1) through 3.1.1(10) in order to learn that B 
is singular. Usually, an examination of the rank and structure of 
3(t^) will suggest the choice of generalized parameters p.^ 

Mesh spacing . - The convergence properties of the TMA are de- 
veloped in Section 8.3. Given the present hypothesis that the gen- 
eralized Jacobian (6) is not singular, it is shown that quadratic 
convergence is obtainable if either the initial guess ip n (t,) is 

U K. 

•JSf 

"close enough" or the time interval is short enough. "Short enough" 
means that the time interval does not greatly exceed the dominant 
local "time constant" of the system. The latter can perhaps be 
judged from the physics of the system, or may be derived from the 
Jacobian matrices 3. 1.1(7), evaluated at particular times of 
interest. 

More practically, a heuristic scheme could be employed to de- 
crease mesh spacing whenever convergence difficulties are noted - 


*C.f. Eqs. 8.3(7) and 8.3(8). 
**C.f. Eqs. 8.3(10) through (12). 
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as evidenced, for example, by an excessive number of trials to reach 
convergence. 

3.1.4 Examp le-Zermelo ' s Problem 

The following example will illustrate the transition matrix al- 
gorithm. It also demonstrates the importance of choosing proper it- 
eration parameters as discussed in Section 8.1.1. As Figure 3-2 in- 
dicates, a boat has a constant speed v relative to the water, 
which in turn is moving with a fixed velocity u in the X-direction. 
It is required to find the heading angle a to minimize the time of 
transit between two fixed points. The system is 


x = u + v cos a y = v sin a (1) 

Upon applying the maximum principle and other steps prescribed above, 
it is readily seen that the optimal control law is 

tan a = (2) 

and the canonical equations are 


where 


X 


U + VTp / P 

y 


vip 2 /p 

♦i 


0 

*2 


0 



(3) 


(4) 


Unfortunately, if ^ and i are taken as iteration parameters, it 
can be shown that the B matrix has a determinant proportional to 




FIGURE 3-2. -ZERMELO'S PROBLEM. 
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~ / p + 1 |~ ^ 1 ^ 2 /P 
" ^2 / p + 1 


( 5 ) 


to the first order in time, so that B is singular and hence the 
pair (ij; , 1 ^) are not proper iteration parameters. On the other hand, 
the pair (a,t f ) results in a tractable result, i.e., to the first 
order 


B = 


_3(x,y) 1 

^ v sin a 

u + v cos a\ 

3(a,t f ) ' 

V v cos a 

v sin a / 


( 6 ) 


whose determinant 

2 

AB = - v uv cos a (7) 

does not vanish if v > u. But, if viu, B is singular at head- 

• r-* 

4 

ings such that 

a = cos 1 (v/u) (8) 

Figure 3-3 illustrates the geometry of the traversal sets % (t^; 
t k’\^ °ear the point (x^, t k >, corresponding to v > u (the oval), 
v = u. (the circle) and v < u (the pointed figure). 

In the latter case, singular B would be experienced if a trial 
mesh point lying outside the envelope were chosen. In that 

case a penalty function formulation as discussed above could be used 
to recover a feasible mesh point. 


3.2 Level II - Iteration in Mesh Point Space, X 


EE 


At this point the selected mesh-points have been connected by 
individually-optimal subarcs, and the machinery exists to re-connect 
them after perturbations. It remains only to perturb the points in 
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3L so that the overall transfer Is optimal also. Clearly, if we 
can show that: (a) the contribution to J arising from each subarc 

is a function of its mesh coordinates only; and (b) define the grad- 
ient of J in I | then all the existing theories of mathematical 
programming may be incorporated in Level II. We demonstrate these 
two crucial points below, before proceeding with the development 

search algorithms in X . 

PP 

3.2.1 Interpretation of the Adjoint Variables as Partial Derivatives 
Theorem . - Consider Problem I, the optimal TPBVP defined on the 
interval <% = [tg,t^], with ? = ?(x,u) only and no P.V.C, terms. 
Assume that : 

(a) The state variable derivatives ?(x,u) are continuous and 
has continuous second derivatives in all arguments and the component 
fg(x,u) is bounded below for all t e 

(b) The admissible controls consist of all bounded, piecewise- 

continuous functions , with finite number and magnitude of discontin- 
uities, whose values u(t) belong to a convex set i) for all 

t e ‘2T; 

(c) The two mesh points x Q and and times t Q and t^ 

defining the terminals of the trajectory are chosen so that discon- 
tinuities in u(t), if present, will coincide with t Q and/or t^ 
Thus, it is sufficient here to consider only single, continuous 
functions with values in fi; and 

* 

This entails no loss of generality since the P.V.C. terms are 
already of the desired form and the neglected parameters g and t 
may be regarded as additional state variables - c.f., 8.1.3. 
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(d) The control u(t) defined over X" and its corresponding 
trajectory x u (t) are unique , feas ible , and optimal . That is, 

—y 

x(t^) = Xq, x(t^) = x^ and the Maximum Principle and the Convex- 
ity, Normality and Jacobi conditions are satisfied for all t £ X 
Then. - (a) the resulting integral criterion value depends upon 


x Q and x^ only, i.e., J = J(Xq,x^) in jKc. S^. is t * ie ^° - 


main of feasible pairs of initial and final points in ^ or 

given tt) ; and (b) the gradient of J(Xq,x-^) exists and is defined 


in R 2N by 


VJCXq,^) = 


- U t Q ) 


+ il>(t f ) 


(i) 


within a scalar multiple. 

Proof . - For condition (a) , optimality implies that an admissi- 
ble control exists in the form U Q ^ t (t) = v(xQ,x^,t), for all t e. ‘Z', 
which solves problem I. Thus, (a) follows immediately upon inserting 
this function into the system equations and the defining integral 
for J. 


The existence of expression (1) is also implied by optimal! ry. 
That is, in the rigorous proof of the Maximum Principle (c.f., 
Chapters 1 and 2 of ref. 137) it is shown that for u(t) and x^(t) 
to be optimal, it is necessary that there exist a non-zero, contin- 
uous vector function *Kt) = (^(t) ,i^(t) , « • • ip (t) ) and a Hamil- 

. _ ^ —y —y * ^-y —y —y 

tonian function jr(4>(t) ,x(t) ,u(t)) = ip ( t ) f(x,u) such that 

and the Hamiltonian is maximized with respect to u e (1 for all t S..2T. 
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Thus, it remains to show that the adjoint variables may also be 
interpreted as partial derivatives of the integral criterion value 
J. This may be readily accomplished, under the present assumptions, 
by simply constructing the first variation, 6J. Since x(t) = 
?(x(t),u(t)) in problem I, we may write 


J = 


f Q (x(t) ,u(t)) + ip T (t){x(t) - f(x(t) ,u(t)) }Jdt 


Recalling the definition of and setting = -1, this may be 
expressed as 


J = 


[-J? + ^ T (t)x(t)]dt 


(3) 


Integrating the last term by parts we obtain 


J = '^ T (t 1 )x(t 1 ) - ^ T (t 0 )x(t 0 ) - 


[ J^+ |i T (t)x(t) ]dt (4) 


•o 


The variation 6J due to perturbations in x(t) and u(t) may now 
be calculated by differentiating under the integral sign in (4) : 

<5J = 4» T (t 1 )6x(t ;L ) - ^(tg) 6x( tg) 


*1 r_ 


^+i) Sx - 6u 


3x 


3u 


dt 


(5) 


0 «- 

(This is justified since ?, and hence have continuous second 

derivatives. Note that if any discontinuities in u(t) are present, 
they occur at mesh points only and do not affect the preceding 
results . ) 
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The first term in the integral vanishes identically because of 

Eq. (2) . And the second term vanishes if the maximum principle is 

3 -jf? 

satisfied. That is, =0 if no constraint is binding and 

6u = 0 if a constraint is binding (otherwise, either the maximum 
principle would be violated or u would be inadmissible). Thus, 

<5J = ^ T (t 1 )6x(t 1 ) - ^ T (t 0 )6x(t 0 ) (6) 

Equation (1) follows immediately from this result. In a similar 
fashion, it may be shown from Eq. (4) that 
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(b) The adjoint vector iKt) also represents the outward normal 




to -si (xQ,tQ;t), the set of attainable states defined in the N + 1 - 
dimensional space [J.x^ • ‘x^ - c.f., reference 137. 

This fact provides a geometrical interpretation of the preced- 
ing results and a means of understanding the role of the second- 
order hypotheses of the theorem. The notation and precise defini- 
tion for the sets of attainability are given in Section 8.1.1. 

These concepts are illustrated in Figure 3-3, where x^ = J is 
plotted against a typical state variable x^t^). The set ^(x^ , t^ ; t^) 
and its boundary are shown in the upper part of the figure. It's 
projection into state space, j^(xQ,tQ;t^) is shown edge-on as the 
heavy line along the x^-axis. The adjoint vector ^(t^) and tan- 
gent plane Tt(t^) are shown at the terminal point x^(t^). By op- 
timality, x(t^) €- 3 ^ and, to the first order, optimal perturba- 

tion lie in it. Clearly, if 3 5& has the smooth, regular struc- 
ture implied by the figure, the first order change 6J due to the 


perturbation 6x^(t^) is given by 

W / 

<5J ; 6x,(t 1 )[or, 6J = + 

1 1 l 


*i <t 0> 




Choosing = -1 as a scale factor and passing to the limit 
6x^ ■+ 0, we recover Eq. (1). 

Evidently this result depends only upon 3^' having a well- 
behaved structure. And, although a rigorous discussion is not in- 
tended, the second-order hypotheses of the theorem actually imply 
proper structure of . Specifically, the continuity properties 



Criterion value Xq or J 



Terminal coordinate value X|(ti) 


FIGURE 3-3 . - GEOMETRY OF THE ATTAINABLE SETS 4 ^ c r n+1 , AND 4 cr n . 
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together with the convexity condition mean that does not have 

"corners" or other irregularities where ^(t) would be undefined. 

The normality condition means that x(t^) ^ (c.f., points "a" 

and "b", at which = 0) ; and hence that the scale factor choice 

= -1 is legitimate. Finally, the Jacobi condition implies that 

8^ is not "folded" or multi-sheeted in a neighborhood of x^ , thus 

establishing uniqueness in a local sense. 

3.2.2 Necessary Conditions in X 

: EE 


For an individual arc the minimum criterion value J, . is 

k,min 

unique if it exists at all, and it depends only on the end-point co- 
ordinates as shown above. Hence, the original integral criterion 
may be expressed as 

K-l 

J ' Ty J k.»in < VVi > (1) 

k=0 } 

Thus, the upper level problem reduces to a problem in ordinary cal- 
culus or mathematical programming, namely: 

minimize JCxqJx^ • • ;x K ) (2) 


where 


If X is the whole space 

pp 


r€l 

It pp 

the classical necessary conditions 



K A 



positive 

definite 


( 3 ) 


apply to the problem thus stated. 
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More generally, the celebrated Kuhn-Tucker conditions, c.f., 
reference 45, would apply if X were a general convex subset of 

pp 

But for convenience, let us first consider the case where the 

points in X are unconstrained. 

PP 

Unconstrained mesh points . - In view of theorem 3.2.1, the 
above first partial derivatives are given by 




-* (t k > 


In we may, therefore, write the gradient vector as 


-V J = 
X 


^(t 0 ) - 0 
^(tp - <Kt*) 


0 - ip(tg) 


At a given point in time, the components of the gradient of J are 
defined in the space of mesh points; hence the name, "mesh gradient", 
for the present concept. Assuming that exists, the Hessian or 

second-partial matrix K may be determined by re-arranging Eq. 3.1 


(3) in the form 


(t k ) 


E 1 F 6x(t k ) 


g|h <Sx(t. . n ) 


E k = - W 


f, = K 1 

k k 


G. = c v “ D i_E. A, 
k k k k k 


«k = °k B k 
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Then by making the following substitutions: (a) 3.1 (6) into (7); 
(b) (7) into (6); and (c) (6) into (5); collecting terms, and for- 
mally differentiating the result, it can be shown that K has the 
block-tridiagonal form 


( 8 ) 

K-l 

L_ 1 U K-1 | "K-l 

Equations (5) and (8) may be used to verify the necessary conditions, 
such as (3), which apply to the upper level. 

P.V»C Terms . - Since J is in general a sum of integral and 
point-value contributions, the general gradient is obtained by add- 

3 p k ( V 

ing a term such as - to each element of the gradient shown 

8x k 

in Eq. (5). That is, if p fc = p^x^) only, then 

= * € *NK <9) 


V J = 

3 





where 
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+ + . 3p o (x o ) 

V -*o + ~ 7 

9x_ 



-*■ -H" 

*k’"V *k + 


3 p k ( v 

3x k 




(10) 


and 


9 P k 

Similarly, a term such as — • must be added to each diagonal block 

9 \ 

in the second derivative matrix of (8) . 

In the more general case where the kr^- PVC is influenced by 
mesh points other than xi^, the additive terms indicated above are 
replaced by summation over the appropriate indices. That is, if k 
represents the range of indices for which p^* is influenced by x, , 
the general gradient term is 


*K ~ 


*K + 


VV 

3x k 




k* 3x k 


(U) 


The second-derivative matrix K is similarly modified, keeping in 
mind that an off-diagonal block will be created each time k* ^ k. 

Constrained mesh points . - A generalization of the transversal- 
ity condition may be derived by noting that any feasible perturba- 
tion 6x fc must be in the hyperplane tt^, tangent to at the 

. T 

point x^. i.e. 



(12) 


for all linearly independent vectors 



in X . 
k 


That is , the 
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gradient components associated with x^ must be collinear with the 
gradient to at x^. More general conditions derived from 

Kuhn-Tucker theory are presented in Appendix B. From henceforth, it 
will be assumed that all these conditions have been enforced, so 
that the gradient is understood to be taken with respect to feasible 
directions only. 

3.2.3 Direction-of-Descent Algorithms 

Having now developed the structure of the criterion value func- 
tion and the necessary conditions that apply to its minimization, it 
is finally appropriate to consider numerical algorithms for implem- 
enting Level II. There are many which could be considered, for in- 
stance references 44, 45, 48, 49, 57-60, 67, 68, 83, 87, 122, 138, 
176, and 185. Three algorithms of present interest are presented in 
Appendix B and their characteristic features are developed. These 
are: 

(a) An apparently - novel, ’'maximally smooth" linear step 
algorithm ; 

(b) A second-order step algorithm; and 

(c) A modification of the Fletcher-Powell algorithm in which 
the unidimensional-search step length is computed in closed 
form from second-order data (as opposed to numerical linear 
searching) 

These algorithms differ significantly in their convergence rates 
and other details, but all function by defining and searching along 




successive directions of descent in X . Hence, it is of interest 
PP 

to consider the convergence properties of direction-of-descent al- 
gorithms as a class. In the present terms it is possible to define 
candidate directions of descent as 

s = {unit vectors in X | s • ¥ < 0} (1) 

That all such directions are directions of descent may be shown as 
follows. Suppose that the NK x NK second partial matrix 


3 JO*) 


3X.3X. V X (X * ) 
i J 


( 2 ) 


is well defined and has a uniform bound in some neighborhood of the 
vector 

L*i "" *kJ =z * 

(By uniform bound, it is meant that there is a number . m > 0, such 
that 


|z T [V 3 .4'a*)]^|< m [ z| 2 (3) 

for all z <£ R^. ) 

Then, use Taylor's theorem with remainder to express the value 
of J resulting from a step of length a > 0 in the s-direction, 
from X*: 

JOT* + as) = J(X*) + us T V x J + ± a 2 s T [V s n£)]s (4) 


where 
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0 <. 0 <. 1 

Now, recalling that s is a unit vector, defining 3 to be the neg- 
ative of the cosine of the angle from s to ¥ and using the uni- 
form bound relation we get 

2 

J<2* + as) < J(X*) + aB | V^J | + (5) 

Since 3 < 0 by hypothesis, it is always possible to find a suffi- 
ciently small a > 0 such that the right hand side of (5) is less 
than J(X-*)(e.g., any value of a in the range 0 < a < -jg- |sVjJ|). 
Hence, the vectors s defined above are directions of descent. 

A direction of finite descent is similarly defined as 

s = {unit vector in x| s • Y = < -e < 0} (6) 

At this point it is intuitively clear that descent algorithms 
based on always searching along a direction of finite descent, will 
eventually converge if J is bounded below. This is of fundamental 
importance and is therefore stated" as a formal: 

Theorem . - Given the hypotheses of section 8.1.1, a uniform 
bound forEq. (2) above, and that J is bounded below. Assume that 
an algorithm using only directions of finite descent is applied to 
the function J (X) . 

Then . - The descent process will eventually converge, in the 
sense that J will ultimately approach its "nearest" local minimum. 
Proof . - From the triangle inequality and the above mentioned 


bound, it is seen that 
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JO* + as) ^ J(X*) + a [j 2 - - | $Vjjj (7) 

where a, g, and m are as used above. 

Arbitrarily set 

« = ~|P VJ I (8) 

then 


Recalling that 


JO*) - J(X* + as) 


> l ^v . J | 2 

2m 


(9) 


3 e (10) 

which follows from the definition of a finite descent direction, 
consider the sequence {.X^} and assume it has a limit point X. Form 
the sums 

P-1 p-1 2 

C |J V - J< W1 - J( V - J «p> >7-1) l«7 J J an 

fco 2m V 


(Eq. (11) follows from the fact that Xj^ = X k + as k ) . Thus, 


p-1 

j v ij< v 

k-0 


'V J k 


( 12 ) 


Now, taking note of the fact that J is bounded below (because f Q 
is bounded below and JT is finite) and passing to the limit 


u- J« P > * JO,) - £ £>k ,J k 
P 0 



(13) 


then shows that the series 


converges, which means that 
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lim |g VJ | 2 = 0 (14) 

i it K 

k^*» 

This implies that 

lim | eVJ | 2 = 0 (15) 

k-*» k 

because of (10) , and hence 

lim | VJ | ^ = 0 since e > 0 (16) 

k-x» 

by definition. Thus, any algorithm which produces directions of de- 
scent will eventually converge. 

In passing, it should be noted that the finite descent condi- 
tion (6) is viewed here as an empirical, verifiable restriction on 
admissible search directions, not necessarily as an analytical prop- 
erty of particular algorithms or objective functions. The class of 
"gradient-restoration" techniques (c.f., ref. 122) functions by, in 
effect, keeping track of |g] and restoring the gradient direction 
(|g| =1) as the current search direction whenever a predetermined 
critical value of e is reached. 

To summarize, the essential Level II operations are to compute 
the criterion value J(x Q , 2 ^- • -x^ the S radient vector (3. 2. 2(5)), 
the Hessian matrix K (3. 2. 2(8)) and then apply a direction-of- 

descent algorithm to find the optimum mesh point locations in X . 

PP 

The overall MG method then entails the following major steps. 

(a) An initial set of mesh points is selected by judgment, ob- 
serving the criteria of 2.4(c) and 3.1.3. 

(b) The resulting sequence of optimal TPBVP's is solved using 
the TMA (summarized at the end of Section 3.1.2), or any other 
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appropriate method. This comprises one cycle of Level I and pro- 
duces the data, needed in Level II for step i. 

(c) Level II then computes the quantities and K^, 

and applies one of several appropriate search algorithms to define 
a direction of descent s. in X and a step length a. such 

i pp * 6 l 


that 


X. , , = X. + a.s. 
i+I i ii 


is minimized with respect to , and 


J..., < J. 

l+l l 


if a, > 0 

l 


(d) The necessary conditions applicable to Level II, c.f., Eq. 
3. 2.2(3) are checked at each step. Steps (b) through (d) are re- 
peated until 3. 2.2(3) is satisfied within acceptable tolerances. 



4. SAMPLE APPLICATIONS 

To illustrate the previous developments, the Mesh Gradient ap- 
proach has been applied to four sample problems. 

Example 1 „ a minimum effort control problem with field-free 
dynamics, is intended merely to . illustrate the Mesh Gradient method, 
using several different descent techniques. 

Example 2 . involves a conjugate-point problem previously ana- 
lyzed in reference 24, and illustrates the present method's strong 
tendency to reject conjugate solutions. 

Example 3 . the optimal control of an unstable Van Der Pol os- 
cillator, is pursued in considerable detail to demonstrate the effi- 
ciency of the method for highly unstable problems with numerous in- 
termediate constraints. Comparisons with alternative methods of 
solution are included. 

Example 4 , illustrates the method's ability to efficiently 
solve optimal, multi— impulse and low— thrust space trajectory problems 


and develops the existence of a hitherto-unsuspected class of 
solutions . 

4 . 1 Example 1 - Minimum-Effort Control 

As a simple illustrative example, consider the following problem 

>10 


Minimize J = — 


u (t)dt 


( 1 ) 


with 


x = u; x(0) = x(10) = 0 
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Obviously, we have that = ip, and ip = 0. Hence, optimal tra- 

jectories consist of straight lines, with slopes given by the local 
value of ip, between successive terminals. The components of the 
vector Aip referred to previously, are given in this problem sim- 
ply by the change in slope at each interior mesh point. 

For the sake of being definite, let there be 11 equally spaced 
mesh points, initially located as follows: 

0 ° o 

x 0 = X 10 = ° 


X i = 


1 < ii9 


( 2 ) 


The corresponding values of are all zero except for 

= Aljjg = - 1. 

4.1.1 Descent Via Programming 

Proceeding by steepest descent we obtain the sequence illustrated 

below. (Note that because of symmetry about t = 5 , we have x^ = 

Xg, etc..) In Figure 4.1, the initial guess x and the corres- 

ponding Aip are shown in heavy lines , the succeeding state- 

iterates are indicated by light lines and identified by step number. 

(Aipj~, but no further Aip k 's is shown.) It is apparent from the 

figure that convergence is very slow. That is, the error measure 

I 

maxlx^Ct) - x 0 pt(t) is not reduced appreciably until 9 steps have 
t , k 

been taken (or in general, K steps if K rather than 9 interior 
mesh points had been used) . s 

This slow convergence is a well-known characteristic of the grad- 
ient method. However, if the conjugate-gradient algorithm is used 



State-variable iterates, Gradient iterates, A^j 
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0 1 2 3 4 5 

Time 


FIGURE 4-1 . - INITIAL-GUESS DATA 
FOR PROBLEM 4. 1. 
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instead, exact convergence is attained after 5 steps (or about K/2 
steps generally), as Table 4-1 shows. 

Table 4-1 

Conjugate Gradient Step Data for Problem 4.1 


Step 

0 

1 

2 

3 

4 

5 

J 

1 

1/2 

1/3 

1/4 

1/5 

0 

X 1 

1 

1/2 

1/3 

1/4 

1/5 

0 

X 2 

1 

1 

2/3 

1/2 

2/5 

0 

X 3 

1 

1 

1 

3/4 

3/5 

0 

X 4 

1 

1 

1 

1 

4/5 

0 

X 5 

1 

1 

1 

1 

1 

0 

S 1 

-1 

-1/4 

-1/9 

-1/16 

-2/25 

0 

S 2 

0 

-2/4 

-2/9 

-2/16 

-4/25 

0 

S 3 

0 

0 

-3/9 

-3/16 

-6/25 

0 

S 4 

0 

0 

0 

-4/16 

-8/25 

0 

S 5 

0 

0 

0 

0 

-10/25 

0 


Symmetrical about x,. , s^. 


The preceding transparently simple example does more than merely 
demonstrate one method of applying state-gradient theory to optimiza- 
tion problems; it also illustrates the primary disadvantage of that 
method. I.E. because of the fixed boundary points, the initial step- 
lengths are relatively small and convergence is correspondingly de- 
layed. It is clear from Figure 4-1 especially that there is no sig- 
nificant reduction of max|x^.(t) - x opt; (t)| until the boundary-point 

b s j 

influences have propogated clear across the trajectory. 
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4.1.2 Descent by First or Second Order Step 

Anticipation of the above difficulties led, in Section 9.3.1 of 
Appendix B, to an attempt at constructing a direction of descent s, 
which vanishes in a smooth, continuous fashion at the fixed termi- 
nals. 

For the same problem and initial iterate considered above, the 
g, h, s, and functions of Section 9.2.1 may be computed as 

follows: 

J = 1 
h = x 


ah 

8x 


0 


ah 

at 


o 



x (t) : as above with 

x* = - 6 (t - 1) - 6 (t - 9) 

g = 0 

Aip = + 6 (t - 1) + 6(t - 9) 
Hence, search directions are defined by 


s = X (6 (t - 1) + 6 (t - 9)) 
s = + X(u(t - 1) + u(t - 9)) + a 
s = + X (r (t - 1) + Q(t - 9)) + at + b 
Here, <$( ), u( ) and r( ) denote the unit impulse, step and ramp 
function of the indicated argument; a and b are constants of inte- 
gration. Clearly b = 0, a = -X, hence 

s = X(-r(t) + r(t - 1) + r (t - 9)) 
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With X = 1, s is precisely the negative of Xq, and the optimum, 
x(t) = 0, is attained in one step. It is apparent that one-step 
convergence will always be attained by this method in the special 
case where the state equations are field-f^ee (the x^ not appear- 
ing at all) and linear in the control. The second-order step method 
will obviously yield equivalent results for this particular example. 
It will clearly provide one-step convergence for a general linear 
state, quadratic criterion system, since the variational equations 
of Section 9.2.1 are then collectively equivalent to the "exact" 
transfer matrix of the linear system. 

2 Example 2 - A Problem with Conjugate Points 

As initially stated, the preceding example serves merely to il- 
lustrate the application of the Mesh Gradient technique with several 
descent techniques. Before proceeding to more complex examples it is 
perhaps of interest to provide a simple demonstration of one of the 
method's more important qualitative features - namely, its ability, 
for some class of problems, to reject conjugate-point solutions. 

A conjugate-point problem which may be solved analytically ap- 
pears in Bryson & Ho (ref. 24); it is: Minimize the time required to 

reach (x^,0) from the origin for the system 

x = V cos 9 

I 

y = V sin 0 

where the scalar velocity is 

V = V Q l/l + y 2 /h 2 

and 0 is the scalar control direction. 
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Upon introducing the variational adjoint function and Hamilton- 
ian, replacing t by x as the independent variable and using the 
maximum principle, we obtain the canonical system 



tan 8 


J 


sec 8 
v 


( 1 ) 


2 2 
(IT + y 2 ) 

with boundary conditions 

y(0) = 0 y(x f ) = 0 (2) 

In deriving equations (1) > use was made of the optimal control law 

\p = - tan 8 [cos 8/v] (3) 


where the bracketed term is constant by virtue of Snell's law. Bry- 
son & Ho showed that this problem can be solved explicitly in terms 
of standard elliptic integrals, and their solution (which did not 
include the heavy "envelope" curves) is presented as Figure 4-2(a). 
This figure shows some of the minimum time trajectories and contours 
of constant transfer time. Clearly there is a conjugate point 

X 

at = tt, y = 0; the trajectory defined by 8 = 0, is optimum for 

x f x f 

0 < — <. tt but not for -r- > it. 

V 

Note, that the contours of constant — - on the sketch develop 

an infinite curvature at the conjugate point (i.e., — Fur- 

ay 2 


thermo re, on 
of constant 


y = 0 beyond the conjugate point > it), the contours 

V ... ' 

k - have a discontinuity in slope upon crossing the 


X-axis . 
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Although this complete solution is already available, it seems 
that some further discussion might be instructive. Note, first, 
that the conjugate point appearing at x^/h = tt is actually the 
"leading edge" of a cusp-shaped region filled with conjugate points. 
The envelopes delimiting this region are shown by the heavy black 
curves in the figure. Points within this region may be reached in 
at least two ways. This may be understood when it is realized that 
the solution curves for each value of 0 q are periodic; so that, a 
given point "P" may be reached either directly or by a trajectory 
containing one or more intermediate zero-crossings. See Figure 

i 

4-2 (b). 

Of perhaps greater interest is the fact that the Mesh Gradient 

approach appears to be unaffected by the presence of these conjugate 

Ax 

points if the mesh spacing is taken to be — < u. This may be un- 
derstood with the aid of Figure 4.3, whep it is recalled from Eq. (3) 
that 

ip * - tan 0 (4) 

Let us assumed, for simplicity, that a single, symmetrically spaced 

mesh point is imposed at = x f /2 . The perturbed trajectory is 

constructed by running arc 1 from 0,0) to (y^,^) and arc 2 from 

(O.x^ back to (y 1 ,x 1 > . The gradient of J at (y 1 ,x 1 > is clearly 

proportional to 2 tan 0^. Thus, J will decrease when y moves into 

the acute angle formed by the 2 arcs. If x f < tt, the acute angle 

is always on the side nearer to the X-axis (Fig. 4-3 (a)) so that 

proper numerical searching will then select y = 0, 0- = 0 as the 

optimal trajectory. 




FIGURE 4-3. - DIRECTIONS OF DESCENT FOR PROBLEM 4. Z 
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f 

On the other hand, if > it, (Fig,, 4-3 (b)) the acute angle is 
initially away from the X-axis, if y^/h is small, Hence, searching 
will proceed in the direction of increasing y^/h until an optimum 
height is found. 

Although the MG technique’s ability to reject conjugate solu- 
tions for some problems is recognized as a valuable feature, it. 
clearly does not hold true for all problems and all perturbations. On 
the surface of a cylinder, for instance, there are two geodesic 
curves joining each distinct pair of points, and in contrast to the 
example above, these cannot be embedded in a single, continuous fam- 
ily of varied trajectories. On the other hand, the Jacobi condition 
may be verified for individual arcs as shown in Section 8 . 2.3 and 
8.2.4. Moreover, the condition can be verified for a total problem 
if the overall transition matrix is known, (It may be found by 
merely forming the product of subarc transition matrices.) The re- 
sulting A,B, C, and D blocks may then be used to verify conditions 
8.2(2) and 8.2(3). 


4.3 Example 3 - Control of an Unstable Van Per Pol Oscillator 

The basic computational steps of the Mesh-Gradient method, out- 
lined in Chapter 3, were applied to the following sample problem. 


where 


Minimize J 



(x l + 


+ u )dt 


( 1 ) 


X x = (1 - x 2 )x 1 - x 2 + u; x 2 = Xj^; 
x^(0) = 0; x^(10) is open; 


x 2 (0) = 3; x 2 (10) is open 
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and u is the scalar, unconstrained control variable., These equa- 
tions represent a nonlinear oscillator which displays limit cycle 
behavior in the absence of control. The reasons for selecting this 
example were that : 

(a) Its low dimensionality facilitated checking, debugging 
and modifications; 

(b) By virtue of its nonlinearity and instability it provides 

a realistic exercise for comparing optimization methods; and 

(c) It allows direct comparison to be made with other techniques 
because it was also considered by Lasdon et al (ref. 100) 
and Mitter (ref. 123) in their studies of control-iteration 
methods. 

4.3.1 Computing Codes 

By introducing adjoint variables and xp^ and defining the 

Hamiltonian function as 




L 

2 \ 


, , / 2 

2 2\ 

1 

- x„ xJ - 

x„ + u 

+ ii»_x, - X, 

+ x_ + u ) 

\ 

2 1' 

2 

y 2 1 VI 

2 / 


(3) 


the steps called for in Chapter 3 may be carried out straight for- 
wardly. That is, the adjoint equations are found by differentiating 
Eq. (3), and the maximum condition yields 


'Pi 


u 


(4) 


opt 2 

which is then used to eliminate the control from both the state and 
adjoint equations. The resulting forms are finally differentiated 
with respect to x^^ and xp ± to obtain the differential equations 
defining the transfer-matrix elements. 
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4. 3. 1.1 Descent by Modified Fleteher-Powell Program 

The first option involves using a modified Fleteher-Powell tech- 
nique to guide the upper-level search and a linearized initial-value 
method for the lower-level boundary problems « The organization of 
the code is indicated by Figure 4-4 (a) and outlined in Section 9„5„ 

A more-or-less self-explanatory program listing is available from the 
author o The Fleteher-Powell routine, "FLPL", was adapted from one 
furnished by Prof. J. D. Schoeffler (CWRU) . The modification refer- 
red to above consisted mainly of using Eq . (9.4(4)) to predict the 
descent step length (by means of subroutine "STEPL") . (The undimen- 
sional numerical search capability of the original routine was re- 
tained as a backup feature, however.) 

4 0 3„1„2 Descent by Second-Order Steps 

Also studied, was a version of the MG method in which a second 
order step was computed as per Section 9.3* The organization of the 
code is illustrated in Figure 4-5, and the program listing is also 
available from the author. Recall that the matrix inversions re- 
quired for lower-level purposes serve double duty here; only one 
more N x N matrix inversion is required to accomplish the recursive 
upper- level solution. 

4.3.1.3 Generating Initial Guesses for Lower -Level Iteration 

Both codes contain a subroutine, "GUESS", whose purpose is to 
predict the change in the initial \p values associated with each 
descent step. This keeps the lower level functioning very 
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FIGURE 4-4. - PROGRAM ORGANIZATION, FLETCHER -POWELL METHOD. 










Print 

output 


FIGURE 4-5. - PROGRAM ORGANIZATION, SECOND ORDER METHOD. 
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efficiently once an initial feasible case has been found. The fol- 
lowing is a systematic way of generating an initial feasible case: 
Differentiate the state equations once with respect to time and 

eliminate appearances of adjoint derivatives (such as by sub- 

* 


stituting the adjoint differential equations • Solve for the 

in terms of the x. , x„, and x.. These values (ip.) are the re- 
j 3 3 i 

quired initial guesses on the adjoint variables. They are such that 


the first-trial variational subarcs are tangent to the initial- 
s'* 


estimate state trajectory at the selected partition points. For 

the present study, the initial feasible trajectory together with its 

slopes and curvatures was determined from the equations 

2 \ 


= 3 exp 


t_ 

10 


o 

<1 


-6t 

10 


("io”) 


Essentially, this equation is arbitrary. It was chosen on the basis 
that it satisfied the boundary condition of (2) and seemed to yield 
a plausible type of asymptotic behavior. 

4.3.2 Numerical Results 

The problem discussed above was programmed for numerical solu- 
tions as shown in Figures 4-4 and 4-5, and was then investigated with 
the aid of a digital computer. A first and general observation result- 
ing from these computations is that the canonical equations tend to be 


In which the maximum condition has been used to eliminate u. 
** 


Here assumed to be represented by the first three terms .of its Tay- 
lor expansion near each mesh point, i.e., x° , x° and x°. 
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very unstable (as previously anticipated) . This was due in part to 
the fact that the fundamental matrix, if evaluated along the entire 
trajectory, contained some extremely large elements (in the "B" 
block particularly). Additionally, nonlinear effects cause these 
same elements to grow with extraordinary speed as the result of per- 
turbations. The two effects combined to make this problem completely 
intractable for initial-value iteration methods when it is treated as 
a single arc. Even when initial values resulting from previous, con- 
verged multiple arc solutions were used, the perturbations propaga- 
ting from the mesh-points (due to numerical boundary-value tolerances 
_6 

on the order of 10 ) invariably cumulated in a machine "overflow" 

and rendered further computation impossible. This general type of a 
nonlinear instability, although perhaps somewhat exaggerated in the 
present example, is often encountered in aerospace trajectory optim- 
ization problems, Initial trials quickly established the fact that 
either of the present Mesh Gradient algorithms will reduce the sta- 
bility problem to insignificance if at least 10 subarcs are used. 
Routine solutions, even from random initial-guess trajectories, can 
then be made. Six subarcs are sufficient if even the slightest ef- 
fort (see Eq. (6)) is made to obtain a plausible-looking initial tra- 
jectory. Using 4 subarcs or less led to persistent difficulty even 
when a very good initial guess was used. Five subarcs seem to repre- 
sent a marginal case with results (i.e., convergence or overflow) de- 
pending upon the fine structure of the initial-guess trajectory. 

These considerations led to the selection of 10 subarcs as standard 
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for the numerical calculations reported below. It was also found 
that at least 30 integration steps were necessary (e„g., 3 per sub- 
arc) to produce 4 digit accuracy in the results. Except where other- 
wise noted, 5 steps per subarc were used herein. 

4. 3. 2.1 Comparison With Previous Methods 

The critical question is whether the present approach repre- 
sents any improvement over the most appropriate existing method. 

Since initial value methods proved quite impractical for the chosen 
problem it appears that some type of control iteration method should 
be the basis of comparison. These avoid the computational stability 
problem by dividing the integration into forward and return passes. 
The forward pass involves state variables only, using a predeter- 
mined control function; this is a stable process. The return pass 
consists of integrating the adjoint equation; this is stable also, 
because the integration is backward in time . The most recent and 
very probably the best methods in the control iteration category are 
the conjugate gradient and second variational techniques discussed 
in the paper by Lasdon et al (ref. 102). Some computational results 
from that paper are duplicated in Figure 4-6 and compared with the 
present Mesh Gradient results . 

Rates of Convergence 

In Figure 4-6 (a), the value of J is plotted against iteration 
number for both the state-gradient results and the previous control 
iteration results, It is clear that both the Fletcher-Power and 2nd 
order versions of the MG approach offer more rapid convergence than 
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FIGURE 4-6(A). - COMPARISON OF CRITERION VALUES. 



Gradient trajectory, d.Mdu or |a|| 
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the control Iteration methods, especially for the first step. The 
difference is very pronounced in comparison with the gradient and 
conjugate gradient methods; it is considerably smaller but probably 
still significant in the case of the second variational control it- 
eration method. That is, if a tolerance of 0.5% on the value of J 
is accepted as defining "convergence" for practical engineering pur- 
poses, the comparison displayed in the following table is obtained: 

Table 4-2 

Number of Steps Required for 0.5% Convergence 


Method 


Number of Steps 


MG (2nd order) 2 

MG (Fletcher-Powell) 3 

Second Variational (Mitter) 4 

Conjugate-Gradient (Lasdon) >20 (say 30) 

Steep Descent 20 


The comparison is further elaborated in Figure 4-6(b), where the 
gradient trajectories for the conjugate gradient and MG/Fletcher- 
Powell methods are presented. Clearly the MG approach has produced 
a better state of convergence after 3 iterations than the conjugate 
gradient approach did after 20. 

Computing Time 

The comparison between methods is less clear-cut when performed 
on the basis of computer time required to reach a given degree of 
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convergence. Such comparisons tend to be very misleading unless 
great care is used to insure that identical numerical-integration 
methods, tolerances, and computing machines are involved in all 
cases. This of course, could not be done for the Figure 4-6 re- 
sults. Rough estimates, however, can be made on the basis that, for 
a given computer, integration algorithm and number of integration 
steps, the time should vary as 


where 


N 

eq 


x N 


ds 


x N 


bv 


N = no. eqs. to be integrated 
N^ s = no. of descent steps 

and 


N^ v = no. of trials to satisfy boundary 
values (the average value is used) 

For the conjugate gradient method, = 4 (2 state and 2 adjoint 

equations only) . For the second variational and MG methods , there 

are also 16 component equations of the fundamental matrix, so that 

N = 20. 
eq 

The values of were previously shown in Table 4-2. The 

value of is unity for all the control iteration methods because 

the terminal conditions are free for this problem. For the MG 
methods, it was observed that averaged 2 to 3 trials during the 

first descent step, but was limited to unity subsequently. (This 
means the 2nd and later steps were computed on the basis of predicted 
values of etc.). Then the applicable value of N, 


decreases 
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as N^ s increases; a value of 1.75 (corresponding to = 2) was 

assumed here. The results of these estimates are presented in 
Table 4-3 below: 

Table 4-3 

Estimated Run Times for Different Methods 

Method Estimated Run Time 

(arbritary units) 

MG (2nd order) (20x2x1.75) 70 

MG (Fletcher-Powell) (20x3x1. 75) 105 

Second Variational (Mitter) (20x4x1) 80 

Conjugate Gradient (4x30x1) 120 

Again the results tend to favor the MG methods, although in this 
case the margin is perhaps not decisive. 

To summarize this section, it should be pointed out that the 
problem discussed above has free terminal conditions . and is there- 
fore unfavorable to the MG method. It lends itself especially well 
to solutions by control iteration methods, however, because the ap- 
propriate final adjoint values (namely zero) are known in advance. 

By contrast, the basic advantage of the MG methods does not apply 
to the open-terminal case because there is no way to reduce the 
dimension of the upper-level search. 

From these results, it may be concluded that, even for an un- 
favorable example, the MG methods (a) have definitely better overall 
convergence rates than control iteration methods, and (b) range from 
slightly better to competitive on the basis of run time. 



72 - 


4. 3, 2. 2 Constrained Problems 

It was originally stated that the MG methods are primarily in- 
tended for heavily constrained problems, so that even better perform- 
ance should be expected in that case. The effect of various types 
and numbers of constraints will be considered next. 

Effect of terminal and intermediate boundary conditions . - When 
imposed upon the preceding problem, these actually result in improved 
MG performance by limiting the dimensions of the upper level search. 
Point boundary values are incorporated especially readily in the 
MG/Fletcher-Powell formulation. The constrained coordinates are sim- 
ply excluded from the search vector and the remaining ones are re- 
numbered. Typical results are shown in Figure 4-7, where run time 
(measured in the same arbitrary but consistent units shown in 
Table 4-3) is plotted against the number of point constraints. For 
the MG method, the run time was 105 units (as in Table 4-3) for fixed 
initial point only, and decreased gradually as more and more con- 
straints were added. Control iteration methods on the other hand 

show a rapid increase because the factor N, increases from 1 to at 

bv 

least 2 (possibly much more) in going to a 2 point boundary value 

problem. It is probable that even further increases would attend 

intermediate boundary values (in general, an extra set of multipliers 

is involved with each added point constraint) and there appears to be. 

no likelihood of a compensatory decrease in N, . These results lead 

as 

to the conclusion that even if the MG method is inferior for open 
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FIGURE 4-7. - EFFECT OF TERMINAL AND INTERMEDIATE BOUNDARY 
CONDITIONS. 0. 5 PERCENT CONVERGENCE. 
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terminal problems, it will surpass the control iteration method even- 
tually, if enough terminal boundary values (not to mention intermedi- 
ate boundary values) are imposed. 

Effect of state-variable point inequality constraints . - Another 
advantage of the programming formulation of the MG method is that 
state-variable inequality constraints can be incorporated point wise 
with no difficulty. A series of runs, similar to that dis- 
cussed above except with equality constraints replaced by consistent 
inequalities, yielded the following result: In no case was the 105- 

unit run time of the original unconstrained case exceeded. 

Note, that the point-wise inequality constraint is sufficient 
for many applications, either as an approximation or because it is 
already an adequate description of what is required. 

4. 3.2. 3 Effect of Poor Starting Iterates 

The initial approximation used so far is quite simple, and is 
evidently no better than that used for the control iteration methods. 
Nevertheless, the possibility exists that both estimates may have 
been fortuitously good. Therefore, the previous initial estimate for 
the MG method was perturbed by adding to each partition point a ran- 
dom point in the interval ± (10- t). Although results varied some- 

what, as might be expected, no case of convergence failure was 


For example, because of the hazard from solar flares, it may be re- 
quired that the point of closest solar approach along an interplan- 
etary trajectory shall not be less than (say) 0.4 Astronomical Units. 
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encountered. Generally, no more than 2 or 3 iterations were required 
to recover a case as good as the original exponential approximation 
yielded. The remaining steps then approach the original pattern 
(c.f.. Fig. 4-6). A typical case is listed below. 

Table 4-4 

Iteration History From Poor Initial Trajectory 
(Typical Example) 

Iteration Number Value of J 

0 108.9 (initial) 

1 48.6 

2 28.1 

3 23.0 

4 21.8 

5 21.5 (0.5% conver- 

gence) 

4. 3,2.4 Effect of Varying the Distribution of Computational Effort 

It may be noted that the distribution of effort between the lev- 
els is essentially arbitrary, provided the subarcs are short enough 
to avoid instabilities. That is, it is possible to set up the prob- 
lem on the basis of a few subarcs each containing many integration 
steps, or vice-versa. Imagine that this distribution is parameter- 
ized by 6 (where OiS < 1), so defined that 8=0 means the en- 
tire burden is carried by the upper level and 8=1 means that the 
lower level carries the entire load. For the present discussion 8 
is arbitrarily defined by the relation 
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N. 

3 = 1 -^ ( 6 ) 

sa 

where 

N = integer number of subarcs 
s a. 

and 


N = integer number of integration steps 
When 6 = 0, a fairly conventional initial-value method is obtained 
as a special case of the MG approach. The special case obtained at 
the other extreme (8 = 1) is basically equivalent to the generalized 
N ewton-Raphson algorithm. The best computational results, however, 
are obtained for intermediate values of 6. This is illustrated in 
Figure 4-8, where run time (normalized so that the minimum value is 
unity) is plotted against 6 for the original example (free term- 
inal problem). A value of N. =36 was used for all cases, the 
steps being distributed in an appropriate manner between the sub- 
arcs. Clearly there is a ve:ry broad minimum of run time around 

6 = 0.14 to 0.20 for this problem (N = 6 to 8) . For low values of 

sa 

6, the bottom level boundary value iterations are seriously hampered 
by the previously described numerical instabilities. For increasing 
8, these boundary value searches improve rapidly and soon reach a 
point where, for all practical purposes, 1-step convergence is ob- 
tained, Further increases of 6 cannot improve the lower-level op- 
eration beyond this point, but does increase the dimension of the 
upper level search. Thus, the run time gradually increases as g 1. 
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0 .2 .4 .6 .8 1.0 

Level distribution parameter, (3 


FIGURE 4-8 . - EFFECT OF VARYING THE DISTRIBUTION OF 
COMPUTATIONAL EFFORT BETWEEN LEVELS. 
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There are two important points to be noted, from this example . 
First, the best results are found for the general case, and not for 
either of the special cases which existed previously. Secondly, the 
importance of choosing the appropriate value of B has been demon- 
strated. In fact, the example suggests the possibility of devising 
an adaptive loop (in effect, a third "level" of the Mesh Gradient 
method) which would vary B in an efficient if not strictly optimal 
fashion, as the problem runs. This is considered to be a valuable 
potential feature of the present approach and warrants further study 
and experimentation. 

4 . 4 Example 4 - Optimal Space Trajectories 

As a final example, the Mesh Gradient technique is applied to 
the following rendezvous problem in space trajectory mechanics (see 
Fig. 4-9). 


with 


Minimize J = ■rr 

zp 



|a(t) | 2 dt + f X/ 
k=0 


( 1 ) 


x(t) = v(t) and v(t) = - m 2 (*)x(t) + a(t) 


( 2 ) 


where 


w 2 = [x 2 (t) + y 2 (t) + z 2 (t)] 3 


( 3 ) 


x, y, and z are the three components of x in R^; t^ 
are fixed; the initial and final boundary conditions are 


and t 


K 


X<t 0 ) = * 0 
x(t K ) - \ 


v(t 0 ) " v 0 

?(t k ) = \ 
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Intermediate 
impulse AV 2 
(one or more) 



^Avj depart Earth orbit 


>-Low thrust propulsion 
phases 


'-A V 3 arrive Mars orbit 


FIGURE 4-9. - TYPICAL MULTI-BURN SPACE RENDEZVOUS TRAJECTORY. 



FIGURE 4-10. - BOUNDARY CONDITIONS FOR HIGH AND LOW 
THRUST SUBARC. 
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(Additional constraints may apply at the intermediate times t^ ... 
t^_^.) Physically, this system approximates a rocket having both 
high and low- thrust propulsion devices. See Figure 4-10. 

The high thrust propulsion phases are represented, as usual, by 
impulsive velocity changes at K + 1 discrete points along the tra- 
jectory. Each of these is considered to be accomplished in a single 

step, such that the actual velocity resulting at the end of one arc, 
->+ 

v (t^) is instantaneously corrected to equal the nominal velocity 
v^ defined for the k~* mesh point. Thus, 

4? k ‘ ' t(c k ) - \ <5a) 

(The initial velocity for the succeeding arc, v (t^), is set equal 
to the nominal velocity v k , i.e., ^ (t^) = v^. The cost is computed 

|Av k | 2 - !?V k ) -5 k | 2 (5b) 

Note, that the total point-value contribution to J thus reflects a 
sum of squares, as opposed to the summation of magnitudes usually 
considered. This choice emphasizes the effect of large individual 
Av’s and is hence considered more appropriate for the case where an 
individual rocket stage is to be used for each high- thrust maneuver. 

The cost of the low thrust propulsion phase is represented by 
the integral term in Eq. (1). The integral involves |a| rather 
than |a|; thus, we are considering a variable rather than constant 
thrust device. This choice was dictated by convenience; however, it 
is known that variable- thrust performance can be well approximated 
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by as little as 2 or 3 optimal constant thrust levels (ref. 46). In 
addition, methods exist for deriving constant- thrust performance from 
variable- thrust results (ref. 77 and 186). 

Finally, the parameter 6 indicates the relative efficiency of 
the low-thrust and high-thrust systems. The trajectory solutions 
will vary from all-impulsive at 6=0 to 100% low variable thrust 
as 6 -*■ 00 . For intermediate values of 6, the solution is required 
to indicate the optimum split between high^ and low-thrust propul- 
sion. 


4.4.1 Subarc Solutions 

A solution by Picard iteration may be obtained by assuming that 

2 

w in Eq. (2) is a known function of time (initially approximated 

by the method of Appendix C) , solving the resulting linear-quadratic 

2 

problems, computing the new w function by Eq. (3), and iterating 
to convergence. 

When the maximum principle is applied, we obtain (for each mem- 
ber of the sequence) the following canonical system. 

-> -* 
x = v 

2 . 

v = - to (t)x + 6y 

t 2 ( 

X = a) (t)y 

->• 

y = - A 

— y y 

where A and y are the adjoint vectors corresponding to x and 
v, respectively. 


» 
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4. 4. 1.1 Boundary Conditions 

As indicated in Figure 4-10, the applicable boundary and trans- 
versality conditions, for fixed end points x^ and x^^, are 

i(t k> - \ m 0 


v (t k> - \ ■ 0 

* (t k+i> - Vi " 0 


(7a) 




rf E |,?| k ■ 0 

' 1 


k+r k =0 


Or, using (5), 


x(t k ) - ^ - 0 
"■ (t k> - \ ’ 0 


x(t. .,) - 


k+r Vi 


= 0 


(7b) 


-J(t k+1 ) + e(^(t k+1 ) - ? k+1 ) . o 


If on the other hand the end point x^^ as well as the nominal ve- 
locity v k+1 is considered to be variable, the applicable boundary 
and transversality conditions are 


<t k ) - \ ’ 


= 0 


v (t k ) - v k . 0 


K 


- r <w + 


3—1 

2 


E 


9 *k+l k=0 


AvE = 0 

k 


K 




& 


3v ( tk+1 ) 


E 

k=0 


IavI 2 = 0 
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Again using Eq. (5) and the chain rule for partials with respect to 
x, these may be written 


A (t k ) + 8 


Av, 


3Av, 




3Av. 


+ Av, 


k+1 


k+1 


.= 0 


J 


V (t k > + 8Av k = 0 


-A (t, + 8<\ Av, 


k+1 


k+1 


3 Av. 


k+1 


(liV, 


'Vi 


+ Av, 


Vl 


> 


= 0 


Lj 


-M <t k+1 ) + 


6A \+1 = ° 


In applying Eqs. (8b), it is important to note that a variation in 
^ may, in general, change the value of |Av| k _ 1 and |Av| k+1 as 
as well as | Av | k » (Recall Eq. 3.2.2(11).) 

The solution of (6) to (8) may be simplified considerably by 
noting that the three (x,y, and z) components of the solution are 
independent and may be computed separately. That is, a 4x4 funda- 
mental matrix is computed and then applied to each of the 3 compon- 
ents,' as opposed to a single operation involving a 12x12 matrix. 

4. 4. 1.2 General-Case and Low Thrust Solutions in Terms of the 
Fundamental Matrix 

Let us write the 4x4 fundamental matrix as 


a ll 

a !2 | 

b ll 

b 12 

a 21 

a 22 

b 21 

b 22 

— - 

— 

— — 

— 

C 11 

C 12 

d ll 

d 12 

C 21 

C 22 

d 21 

d 22 


( 9 ) 
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where 


2 

-w 0 




0 B 

---- 

0 a)^ 

-1 0 . 


so that 


“VV " 1 




Combining (7), (9), and (11) we see that 


0a 21 C 21 &a 22 C 22 eb 2l“ d 21 eb 22 _d 22 



Or, after solving analytically for 


x(t k> - \ 


(13) 
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and 


v (t k ) - v k 

the unknown adjoint variables at the beginning of the arc are given 
by 


'11 


b 12- 1/S 


L Bb 21 d 21 | 6b 22 d 22 S 


Inverting (14) we have 


Vl- a ll*k' a 12 v k 

'Vl-NV^k - (Sa 22 +c 22 )v kj 

(14) 


X 


— 

_ l 

y _ 

A 


eb 22 _d 22 _1/e 1/e_b i2 


-6b 01 d, 


11 


\+l~ a ll\- a 12 V k 
0V k+l _(ea 21 +C 21 )x k -(ea 22 +c 22 )v k 


(15) 


where the determinant is 

A = b ll (eb 22 ‘ d 22 " 1/3) + (b l2 " !/e) (-6b on + d ni ) 


21 21 " 


When 3 ->■ 00 , it can be shown that X (t ) and V (t ) reduce to 
those values that result from the standard, low- variable- thrust ren- 
dezvous problem> i.e., 


'(O = 


-b 


11 


k b n b 22 b 12 b 21 k+1 21 ^ a 22 V k^ 


21 


+ b ll b 22 “ b 12 b 21 (Xk+1 ' 3;LlXk " 3;L2Vk) 


22 


A (tfc) b ll b 22 " b 12 b 21 (Xk+1 ~ ailXk ” 


d 12 

b ll b 22 " b 12 b 2 1 Vk+1 " a2lXk " a 22 V k ) (16) 
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4.4.1. 3 Solutions for the Case of Impulsive Thrust 

When 8 -*■ 0 in the impulsive- thrust limit, both A(t) and 
y(t) vanish over the entire interval, and no low thrust propulsion 


is used. In this case the initial and final velocities, v (t^) and 

v (t^^) must be chosen so that the desired mesh coordinates 
— y- 

and are achieved at the end points of the trajectory. That 

is, when the applicable upper-left block of $ in Eq. (9) is re- 


labeled as 


a b 
_c dj * 


we have 


a b x 


c d v 


from which we find that the required initial velocity at t, is 

tc 

v (t k ) = b~ 1 (x k+1 - a^) (18) 

and the corresponding (dependent) final velocity at t^ + ^ is 

v 4 " (tj^) = (c - db -1 a)x k - b _1 ax k+1 (19) 

4.4.2 Mesh-Point Iteration 

When the preceding results are substituted into the general equa- 

th 

tion for the gradient vector at a general (k ) mesh point, we obtain 


v- 


A + (t k ) - X (t k ) + 6 


j=k-l 


3Av. | 


3 L 9 \ 


" ( V - 11 <V + s l_j 4v j 
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for the general, high + low thrust case. Notice that a variation of 
a mesh point coordinate will be general change the values of 

|Av|^_ 1 and l^ v l^+l as we -^ as |Av|^. But since a velocity- 


discontinuity is permitted, a change in the desired velocity v at 

K 

y 2 y 2 

a mesh point will affect the values of |Av| k and |Av( k+ ^, but not 
For the low thrust system, the above reduces to 


v- 


X + (t k ) - A-(t k ) 

f(t k ) - T(t k ) 



AX 


Ay 


For the impulsive case, the mesh-point initial velocities are no 
longer independent but must be chosen so that the prescribed coor- 
dinates are attained at the end points of each subarc. Thus, the 
gradient is with respect to coordinates only, i.e.. 


V 


Ui 3 L 3 \ 


As before, the gradient in may be constructed by simply. ad- 

joining the column vectors (20, (21), or (22), taking care that the 
summation index j remains in the same range (0 i j i K) as does 
the mesh point index k. 

The required partial derivatives may be readily computed from 
Eq. (19) to be 

k+1 + + 

3Av. r — v 3v . 3 vj 

( 23 ) 


*k jfe-1 [_ 9x k 3x k 
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where 


3Av. 


k-1 


3v 


)x k 


k-1 -1 

= — D 


3x k 


k-1 


and 


“Vl 8 ^+l . ,,-L. 

— T = — “ — = (c - db a). 


5x k 3x k 


3Av k 9v+ 3v“ 

~T“ = — " = db 

3x. 3x, 3 


3x k 3 \ 


+ b 


k-1 


Thus , in X we have 
PP 


TB 1 A) () Av 0 + (C - DB 1 A) Q Av 1 + 0 


V J = 
X 


- B ^0 + 


(DB 1 ) Q + (B ^^A^ + (C - DB 1 A) 1 Av 2 


-VW 2 + I®*" V 2 + (B ‘ 1a) K-1 


‘Vi + (c - db ‘ 1 a) k-A 


+ + < b ' 1a >k 


A^ k + 0 


(24) 


The partial derivatives for the general case, Eq. (20), may be com- 
puted from the relations 
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3Av, 


3v. 


3v, 


- I 


3Av. 


k+1 


9v, 


’Vi 

3v,_ 


3A \ . 8v k 
d \ s \ 


> 


3Av. 


k+1 


3v 


k+1 


5x k 


)x k 


(25) 


where 


3v, 


3v, 


= b 


21 


4- 

3A 


k-1 3v k 


^i + b 

+ 22 


3y 


k-1 


k-1 3 \ 


-N- 

3v k 


Jx k 


= b 


21 


-4 

3X" 


k-i dx k 


1 k-1 

--- 1 + b 

-> 22 


3y 


'k-1 


++ 

*Vi . 

b 21 


3v 


3 \ 

‘+ 


-4 

ax“ 


+ + b 22 


k+1 


= b 


3v, 


21 


k 3x k 

3X k 

— - + b 

»-*■ 22 

k 3 \ 


k-1 3 \ 

-4 

3y" 




k 3 \ 


3y 


4- 

k 3v k 


(26) 


The eight partials of y and A in Eq. (26) are in turn computed 
by differentiating Eq. (15) and re-labeling indices. For example. 
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-+ 



^ gb 22 d 22 1 /3)a 11 

b ll (eb 22 " d 22 " 1/e) + (b 12 " 1 /^(-3b 21 + d 21 ) 


+ 


( 1 /e - b 12 )(6a 21 + c 21 ) 

b ll (3b 22 ~ d 22 “ 1/B) + (b 12 " 1 /^H-6b 21 + d 21 ) 


(27) 


and the other seven are similarly defined. 

The gradients computed as above are finally used to drive a 
Fletcher-Powell minimization routine, c.f.. Appendix B. 

4.4.3 Numerical Results 

4.4. 3.1 Accuracy of the Approximate Solution 

As point out in Section 4.4.1, solutions may be computed by the 
method of Picard iteration. On the other hand. Appendix C shows an 
approximate analytical form for the w function which is "exact" at 
two points of the trajectory and permits a closed form solution for 
the state and adjoint variables. By inspecting the defining Eqs. 
10(4) and 10(5) for to, it may be seen that to never changes sign 
and to itself is monotonic. By contrast, the "real" to function 
(Eq. 4.4(3)) has to = 0 at every apse or turning-point of the tra- 
jectory and cannot be well-approximated by a monotonic function. On 
the other hand, the "real" function can evidently be approximated as 
well as we please by adjoining a number of monotonic segments, each 
having its own appropriately chosen value of K in Eq. 10(5). This 
is illustrated in Figure 4-11 where to functions are plotted against 
heliocentric travel angle for an Earth- to Mars-transf er trajec- 
tory such as the one illustrated previously in Figure 4-9. The solid 



Schuler frequency, 
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u) = (radius) - -^ 

cj per appendix C 



FIGURE 4-11. - SCHULER FREQUENCIES 
(w) ALONG 300°, 250-DAY EARTH -TO 
MARS ORBIT RENDEZVOUS TRAJEC- 
TORY. 
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curve is the "real" function, the dashed curve shows the single seg- 
ment approximation - which as expected is extremely poor - and the 
dotted curve shows a five-segment approximation closely approaching 
the "real" result. 

Clearly, this approach is compatible with the Mesh Gradient 
technique, with m being approximated by Eq. 10(5) on each of the 
Level I subarcs . The descent process is illustrated in Figure 4-12, 
for a low-variable-thrust vehicle flying the same Earth-to-Mars 
orbit rendezvous mission considered above. The current criterion 
value, J n is plotted against the iteration number n, and compared 
with the "exact", numerically integrated result obtained from refer- 
ence 112. The case n = 0 represents a first or preliminary guess 
using only one trajectory arc; as could be expected, it is a poor ap- 
proximation. Nevertheless, when coordinates and velocities from the 
single arc approximation were used to define mesh points for a 5-arc 
approximation, the much-improved value shown at n = 1 resulted. A 
Fletcher-Powell search (initialized by the n = 1 trajectory) then 

I 

resulted in the descent sequence denoted by n = 2, 3, 4, 5, . . . ; with 
convergence obtained, for all practical purposes, by the fifth step. 

The effect of K, the number of mesh-points used, is illustrated 
in Figure 4-13. Here the value of J obtained at the end of the de- 
scent process is shown as a function of K (or the number of subarcs, 
= K-l). Clearly, the original error is drastically reduced by using 
even 2 subarcs, and 5 subarcs yield a rather accurate approximation. 
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Iteration number, n 


FIGURE 4-12. - CRITERION VALUE VS ITERATION NUMBER FOR 
LOW-VARIABLE THRUST, 300°, 250 DAY EARTH -AAA RS ORBIT 
RENDEZVOUS TRAJECTORY. 
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FIGURE 4-13 . - EFFECT OF NUMBER OF MESH 
POINTS ON THE CRITERION VALUE OB- 
TAINED. LOW-VARIABLE THRUST, 300°, 
250-DAY EARTH -MARS ORBIT RENDEZVOUS 
TRAJECTORY. 
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This illustrates the Mesh-Gradient method’s ability to extend the 
validity of a useful but limited-range approximate solution. In the 
present example, this results in an estimated saving in computer 
time on the order of 10/1 compared to "exact" numerical- integration 
methods such as reference 112. (This is based upon a comparison of 
the number of computational steps required for each method; the ref- 
erence computer program is no longer available for direct compari- 
sons . ) 

4.4. 3.2 Low Variable Thrust Solutions 

The above described procedure, using 5 subarcs, was applied to 
a range of Earth-orbit to Mars-orbit rendezvous trajectories. Re- 
sults are shown in Figure 4-14 with J plotted against heliocentric 
travel angle 0 for a travel time of 250 days. The dotted portions 
of the curves denote cases where the vehicle makes less than one 
revolution about the sun; these are referred to as "direct" trajec- 
tories and are of course repeated every 360 degrees. The solid 
curve, however, represents "indirect" trajectories which wind fully 
around the sun before proceeding to the designated rendezvous point; 
these require more than one revolution. 

Such "indirect", multi-revolution trajectories have proven 
troublesome in the past (computationally) due to numerical instabil- 
ities. This is because the multi-revolution trajectories generally 

2 

involve close perihelion passages, so that the true w function is 
very large, for a time. This in turn causes the systems "effective" 
time constant to be short enough compared to typical travel times. 
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FIGURE 4-14. - COMPARISON OF DIRECT AND INDIRECT 
TRAJECTORIES, LOW VARIABLE THRUST, 250-DAY 
EARTH -MARS ORBIT RENDEZVOUS TRAJECTORIES. 
(CIRCULAR, COPLANAR PLANET ORBIT. ) 
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that the results are extremely sensitive to small errors in initial 
estimates. For example, even the 300° example trajectory discussed 
above (noted on Fig. 4-14 by an asterisk) has a minimum solar dis- 
tance of about 0.1 AU (Astronomical Unit), and stays below 0.25 AU 
for about 25 days out of the 250 day trip. During that time, the 
local "time constant" varies from 3 days to as little as 12 hours, 
and small perturbations may become grossly magnified. 

Therefore, many previous trajectory computation methods have 
experienced significant difficulties in solving the multi-revolution 

TPBVP. Reference 112, for example, did not attempt to compute 
trajectory data beyond 330° even though it was intended as a major 
low-variable thrust trajectory data compilation. 

4. 4. 3. 3 Impulsive Thrust Solutions 

At the opposite extreme, we may consider the case where only 

A i -, 2 

the point-value contributions, \ |Av|^ are significant. Then, as 

0 


pointed out in Section 4. 4. 1.1, the low-thrust system is not used at 

all and Eqs. 4.4(6) may be written simply as 

-*• 2 -> 
x = - 0 ) x 

(28) 

-v 2-y 
X = u> X 

The appropriate fundamental matrix for Eq. (28) is also given in 
Appendix C, and the elements of the gradient vector were defined in 
Section 4.4.2. The resulting algorithm was then applied, as in the 
low variable thrust case, to a range of interplanetary transfers. 
Five subarcs were used for the sake of accuracy. Thus, up to six 
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impulses may be used. Descent and convergence properties were found 
to be substantially similar to those illustrated above in Figures 
4-12 and 4-13, and will not be further discussed. 

Figure 4-15 shows results obtained for a family of 250 day Earth 
orbit to Mars orbit rendezvous trajectories. The criterion is shown 
as a function of travel angle as before; note, however, that 

^|Av| k rather than ~ Av | ^ is plotted. This is merely to 

0 0 

facilitate comparison with previous studies, such as reference 178; 
up to 0 = 360° (the largest value considered in ref. 178) the re- 
sults are indistinguishable. Both reference 178 and other available 
methods (e.g., Refs. 63-65, 73, and 74) could in principle be ex- 
tended to operate beyond 0 = 360°, but for practical reasons this 
has not been done and results comparable to the present ones are not 
available. It would appear, however, that the MG approach may have 
some significant computing-time advantage over the previous methods 
since (even for equal rates of convergence) the w-approximate sub- 
arc solutions require a smaller amount of calculations. In any case 
its time advantage over the reference 178 technique, the only one 
for which direct comparison are available, varies from 5/1 to 25/1 
for computing a single trajectory. These numbers, which basically 
reflect the number of iterations required by the previous approaches 
to solve the Kepler time equation, are probably representative of 
newer methods (e.g., ref. 74) as well as ref. 178. 
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FIGURE 4-15. - COMPARISON OF TWO- AND 
MULTI-IMPULSE TRAJECTORIES. 250-DAY, 
EARTH-MARS ORBIT RENDEZVOUS. 
(CIRCULAR, COPLANAR PLANET ORBITS.) 
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The present approach has another, less obvious advantage com- 
pared even to the most recent available method (ref. 74) in that it 
appears to be more capable, if not infallible, in determining the 
optimum number of impulses. This is due primarily to differences in 
the way that "initial guess" trajectories are prepared. In the pres- 
ent approach the analyst is required to sketch or visualize a 
plausible- looking trajectory and select the location and timing of 6 
(or more) impulses. The descent process then eventually rejects un- 
needed impulses by driving them to zero magnitude. In the examples 
shown on Figure 4-15, this led to 2 and 3 impulse trajectories, but 
(for coplanar planet orbits) never 4, 5, or 6. 

Other approaches, by contrast, typically derive an initial 
guess from a non-optimum (but readily calculable) two-impulse solu- 
tion. The "primer test" as developed in reference 65 from the basic 
theory of Lawden (ref. 107) is then applied to determine whether the 

A 

final trajectory could be improved by adding an intermediate im- 
pulse. If so, a third impulse is added and its location optimized 
as in the present approach, e.g., by the Fletcher-Powell technique. 
(The required fundamental matrices may be computed in closed form, 
c • f • » refs. 33 and 164.) Typically, this results in a locally 

A 

Briefly, Lawden^has shown in Chapter 5 of reference 107 that the 
primer vector (y in the notation of this section) is aligned with 
the velocity ai^d has unit magnitude at the location of an optimal 
impulse. If |y| > 1 at some point, this implies that the trajec- 
tory could be improved by adding an impulse at some interior point, 
e.g., the point where |y| maximizes. 
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optimal 3 impulse solution — i.e., whose two subarcs are each locally 
optimal as judged by the primer test. There are cases, however, in 
which 4 or more impulses are known to be optimal (c.f., ref. 39); 
within the context of Earth to Mars orbit rendezvous trajectories 
this may occur when there is a significant out-of-plane motion. 

The cases shown in Figure 4-15 were re— computed for 270° <. 0 <. 
360°, under the assumption that Mars is located 0.5 AU above the 
ecliptic plane. The results, not illustrated, show that: (a) there 

are in general two distinct locally optimal 3— impulse trajectories, 
only one of which would be obtained by reference 24 or any other ap- 
proach relying exclusively as the primer test; and (b) in many cases 
an (apparently unique) 4-impulse trajectory exists which is somewhat 
better than either of the 3- impulse solutions. The physical logic 
for the 4- impulse solution is that energy can be changed most effici- 
ently in a region of high path velocity - i.e.» near perihelion as 
in reference 178 - whereas the plane-change maneuver is accomplished 
more efficiently in a low velocity region - i.e., farther from the 
sun as in reference 47. For this reason, earlier results involving 
3-impulse, broken-plane transfers (c.f., ref. 35) may be regarded as 
suspect. 


Finally, the MG algorithm was applied to a survey of fast 3- 


impulse Earth-Neptune trajectories. 


4-16, are presented in terms of 


xA 

u 


0 


The results, shown in Figure 
Av| versus travel time and 



Criterion value, I^Av^, m/sec 
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FIGURE 4-16. - CHARACTERISTICS OF OP- 
TIMAL 3-IMPULSE EARTH -NEPTUNE 
ORBIT RENDEZVOUS TRAJECTORIES. 
(CIRCULAR, COPLANAR PLANET ORBITS.) 
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angle as before and required under 2 minutes of IBM 709 A computer 
time. One entirely unexpected result from this survey are the "sec- 
ond minima" visible in the upper right (at long angles and short 
times). These correspond to trajectories which nearly graze the sun 
when the interior impulse is applied. Thus, the velocity is extrem- 
ely high at the impulse point and a very large energy increment can 
be gained from a moderate lAv^. The point of this example is that: 
(a) the MG approach is efficient enough to allow the wholesale com- 
putational study of optimal multi-impulse trajectories, on a scale 
previously approached only for 2-impplse non-optimal trajectories; 
and (b) unexpected, hitherto unsuspected results and conclusions 
may sometimes emerge from such a study. 



5. CONCLUDING REMARKS 

Although the above-described analysis and examples are prelim- 
inary in several respects, they support several basic observations. 

(1) First, the MG approach is computationally feasible for its 
intended class of problems. This further implies that the theoret- 
ical structure, developed up to this point, is well founded. 

(2) Secondly, the method has realized its basic objective of 
handling, in a routine manner, problems which are very unstable, 
and/or involve numerous intermediate boundary conditions. Note 
that it is not necessary to guess initial values of the adjoint var- 
iables. 

(3) Also, pointwise state-variable inequality constraints can 
be incorporated in a routine manner with no sacrifice in run time. 

(4) The method has very strong convergence properties in gener- 
al, and the initial rate of convergence is especially remarkable (re- 
call Table 4-4). Ultimate convergence is quadratic, which is at 
least as good as that provided by any competitive methods. 

(5) As to computational efficiency, the MG method is at least 
competitive with existing ones for open terminal problems and evi- 
dently superior if there is an appreciable number of terminal or in- 
termediate boundary conditions. 
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(6) The method is characterized by a flexible level structure 
which can be varied to suit the requirements of a given problem. 

This has a significant bearing on computational efficiency and sug- 
gests the desirability of adding a third, "problem-adaptive" level 
to the MG method. 

Considerable additional work is indicated to confirm and extend 
these conclusions by applying the MG method to realistic examples 
drawn from practice. It must also be recognized that the method is 
presently in an early stage of development. Major generalizations 
and refinements that could be considered include (but are not lim- 
ited to) the following: 

(1) It would be valuable to have a more general procedure for 
generating an initial feasible trajectory when control-variable in- 
equality constraints are present. That is, aside from resorting to 
penalty functions, how do we correct an unfeasible initial trajec- 
tory? The notion of generalized iteration parameters (c.f.. Sec- 
tion 3.1.3) could usefully be further developed - perhaps by extend- 
ing the work of reference 115. 

(2) Further work in the area of state-variable inequality con- 
straints is required to either develop analytical criteria for sat- 
isfactory pointwise approximation to a regional constraint, and/or to 
incorporate the alternative necessary conditions into the existing 
Level I and II structure. The technique of reference 158 for the 
separate computation of arcs lying on a constraint boundary, in 
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particular, appears to be compatible with the present methods. The 
constrained arcs (or families thereof) could be treated as "active 
manifolds" (see 4(a) below). 

(3) Criteria and procedures for implementing the problem- 
adaptive level XII should be developed. For instance, a heuristic 
method might be devised to sub-divide an arc that is "too difficult" 
(as measured by the number of TPBVP iterations required) and also to 
combine adjacent "easy" arcs. 

(4) Extensions of the MG approach to more general problems than 
those illustrated would be desirable. For example: 

(a) Problems with active constraint manifolds (e.g., which can 
produce a discontinuous change in some of the state 
variables) ; 

(b) Branched and/or segmented trajectories; 

(c) Problems involving extremely large dimension, e.g., N = 50. 

(5) There are also numerous detail refinements which can signi- 
ficantly streamline the computations. For one example, it appears 
to be unnecessary to recompute the transfer matrix at every step. 

For another the K-matrix data could be used, after the first feas- 
ible Level I trajectory has been constructed, to predict initial x^ 
values corresponding to a chosen state perturbation Sx^. Based on 
limited experience with the system of Example 4.3, this appears to 
yield a 2 or 3/1 time savings for all Level I steps after the first. 

(6) Based on the encouraging results of Example 4-4, it would 
be worthwhile to further develop the capabilities of the MG method 
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for approximation purposes. For example, it often happens that 
simple feedback algorithms or other approximations can be derived 
which, though perhaps suboptimal, are quite efficient if applied to 
a short enough problem (the "velocity to gain" rocket steering al- 
gorithm is a case in point, the approach of Example 4-4 is another). 
One of these, used in place of the present lower level iteration 
procedures, would greatly simplify the calculations by eliminating 
the need to integrate and store matrix elements. This may be ex- 
tremely significant for high dimensional problems because the number 

of equations to be integrated would then vary as 2N rather than 
2 

4N . It is tempting to speculate that suitable "steering laws" for 
interesting problems could be derived under the conditions used here- 
in. These, combined with the present results could open a completely 
new line of approach to high-dimensional problems. 

(7) While the MG approach was designed primarily with computa- 
tional applications in view, it appears to have at least some theo- 
retical utility which could be exploited. For example, necessary 
conditions for branched or segmented trajectories appear as nearly- 
obvious corollaries of the basic MG results, but would require lab- 
orious derivations if approached from the usual calculus-of- 
variations viewpoint. 
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7. MAIN SYMBOLS 


Matrices 
A, B, C, D 

E, F, G, H 
K 

V V \ 

V w i 
$ 

Vectors 

1 


g 

-> 

k 


-> 

s 

X 


y 

u 

-> 

V 


V 


T 

PP 


N x N blocks of the canonical (state plus adjoint) 
equations' fundamental matrix 

derived N x N matrices 

KN x KN - dimensional Hessian matrix 

auxiliary matrices used in the "H-Process" for 
matrix inversion 

auxiliary matrices used in Varga's recursion 
formula 

2N x 2N - dimensional fundamental matrix of the 
canonical (state plus adjoint) equations 

time derivative of N-dimensional state vector, 
c. f . , Section 2.1 

ditto adjoint vector, c.f., Section 3.1.1 
auxiliary vector used in recursion formula 
search direction in or 

N-dimensional state vector, c.f., Section 2.1 
N-dimensional state perturbation 
M-dimensional control vector 
velocity vector, c.f.. Section 4.4 
control synthesis 

NK- dimensional state vector in the space of mesh 
points 
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Y 

? 

¥ 


V V v \ 


"k 

Scalars 

a 

6 


x, x t , Uy 


m 



NK-dimensional state perturbation vector in the 
space of mesh points 

tangent vector to a surface in 

N-dimensional adjoint vector 

adjoint discontinuity at a mesh point 

NK-dimensional adjoint discontinuity vector in mesh 
point space 

adjoint perturbation vector \ 

auxiliary vectors 
terminal error vector 

L-dimensional vector of design parameters 
auxiliary vector 
auxiliary vector 

(dimensionless) step length in unidimensional search 
auxiliary function or bound 
Lagrangian multipliers 
auxiliary variable 

number of function evaluations per step 
A bound 
A bound 

independent variable (time) 

Hamiltonian function 

functions describing equality constraints 
criterion value to be minimized 
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Subscripts 


i, j, k, 1, m, n 
o 

opt 

f 

u 

des 

act 

ext 

max 

min 

L, M, N 


general indices 
original or initial 
optimal 
final 

due to the control u(t) with t C 

desired 

actual 

extremal 

maximum 

minimum 

dimensiop of parameter, control, or state vectors, 
respectively 


Point Sets 

^ ( V fc o ;t) 

$(t;t f ,x f ) 

F 

sa 

X 

ft 

r 

x 

pp 

*N 


N + 1 dimensional set of attainable states 
c.f. , Section 8.1.1 

reversed set of attainable states, c,f.. Section 

8 . 1.1 

ttl 

k constraint manifold 

the set of feasible subarc solutions 

basic time interval upon which the problem is de- 
fined, i.e. , (t Q ,t f ) 

set of admissible controls 

set of forbidden states 

set of permissible mesh points 

the real N-tuples 

tangent plane to ¥ 


in 


®N+1 
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Other Notations 


( A ) 

(’) 

( ><!> 


vector 
unit vector 
d( )/dt 
d j ( ) /dt J 

transpose of a vector or matrix 


3 


boundary of a point set 


6, A perturbation symbols 

V , V+ indicates gradient operation with respect to the 

variable z or the components of the vector z 



8. APPENDIX A 


DISCUSSION OF LEVEL I - THE TWO POINT BOUNDARY VALUE PROBLEM 
The general problem described in Eqs . 2.1(1) - 2.1(4) may be 
reduced to a form suitable for numerical solution by defining the 
adjoint or costate variables, the Hamiltonian function and then us- 
ing the Maximum Principle to determine the optimal control law in 
terms of state and costate variables. The resulting, 2N-dimensional 
state-adjoint system has mixed or "two-point" boundary conditions - 
half of which apply at the initial time and half at the final. Such 
systems cannot, in general, be numerically integrated in one pass - 
and as pointed out above the required iteration processes are often 
beset by major difficulties. Thus, to support the development of 
the "Mesh Gradient" approach in its entirety, it is appropriate to 
review the theoretical and computational basis of the optimal con- 
trol TPBVP . 

8.1 First-Order Necessary Conditions 

In this section, several standard definitions and theorems are 
presented, without proof,, in a form appropriate to the present work 
and in the interest of completeness and uniformity. 
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8.1.1 Basic Definitions and Hypotheses 

Continuity and differentiability . - The state-variable deriva- 
txve function f(x,u,3,t) are assumed to be continuous, at least 
twice-dif ferentiable in the components of x and u and at least 
once differentiable in t and the components of p. This implies 
that the Hamiltonian function, c.f.. Section 8.1.2, has the same 
properties. Also, the component fQ(x,u,$,t) which defines the 
integral contribution to the criterion value is regarded as being 
finite and bounded below. 

Admissibility and convexity . - The class of admissible con- 
trols, denoted by ft, is the collection of all bounded piecewise 
continuous controls u(t) , with t £ 7, whose values lie in the 
convex set ft eR^. 

Attainable states . - The set of attainable states in R^ is 
defined to be the collection of end points of trajectories 
emanating from the initial point x^Ct^), corresponding to all pos- 
sible admissible controls u(t) e ft. I.e., 

j^(x 0 ,t Q ;t) = (x u (t)|u(s) «. ft, s SC"} ( 1 ) 


The corresponding entity in is formed by treating the integral 

criterion value as an additional state coordinate, say Xq = J. It 
is denoted by the prime notation, i.e., . 

In a similar fashion the reversed attainability set, 

#(t;t k ,x k ) = {x u (t)|x k e jtf(x(t) ,t;t k )} (2) 
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is defined as the set of all points x(t) from which the point x^ 
x(t^) can be reached. 

The notation ) and ) will denote the attain- 

ability sets corresponding to optimal controls . 

It should be noted that the sets j^Cx^t^t) and $(t;t ^ 

are defined at one instant only (namely, t) . It will also be con- 
venient to define "sets of traversal" as the union, over some sub- 
set of , X' y of the sets of attainability. I.e., let 


* •-> 

f ( x o»t 0 ;t ) 


U 


and 


sc[t Q ,t] 


*f (x 0’V s) 




(3) 




sc[t,t k J 


In passing, it may be observed that although & and 
necessarily convex, they have a nesting property, i.e.. 


are not 


s? (Vw 


^ x 0 ,t 0’ t 2^ 


if 


> 


fc l ^ c ^ 


(4) 


with a similar relation holding for the reversed sets . 

Controllability . - The system x = |(x,S,#,t) Is said to be com- 
pletely controllable in the neighborhood of a point x^ e if the 

„ ■k 

traversal set (t;t k ,x k ) of all points steerable to x, (t, ) in 
a finite time interval (t fc - t) contains an open , N-dimensional neigh- 
borhood of x^o It will be termed partially controllable if 
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includes at least one M-dimensional subset Q** which, in turn. 


contains an open neighborhood of 

ted by the following sketches. 

Completely 
controllable 



It is assumed that all systems 
tially controllable. 


x^. These concepts are illustra- 


Partially controllable 



idered here will be at least par- 


Optimality . - An admissible control u(t) is called optimal if 
the statement v(t) ^ u(t) for some finite subset of VC implies 
that J(v) < J(u). For practical purposes, optimality is considered 


to exist when 


(a) The maximum principle and transversality conditions (c.f.. 
Sections 8.1.2 and 8.1.3), 

(b) The convexity or strengthened Legendre-Clebsch condition 
(c.f„. Section 8.2.1), 

(c) The normality condition (Section 8.2.2), and 

(d) The Jacobi or no-conjugate point condition (Section 8.2,3) 
are satisfied. 

The "sufficiency" of these conditions is treated in references 4, 24, 
and 108, for instance. 
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8.1.2 The Maximum Principle for the Basic Optimal Control Problem 
Consider first a simplified version of the problem defined by 
Eqs . 2.1(1) through 2.1(4), with autonomous dynamics, boundary con- 
ditions consisting of two fixed end points and integral criterion 
only (no PVC) . I.e., Problem I : 


minimize J = / f Q (x(t) ,u(t))dt 


( 1 ) 


with 


subject to 


with 


and 


x(t) = l(x(t) ,u(t) ) 


u ( t) £ fi CRy 


x( V = x 0 e ®N 


x(t f ) = x f 




( 2 ) 


(3) 


(4) 


Let fg>f-L ’ * * * f N be continuous functions on the interval (t^jt^), 
also possessing continuous first partial derivatives with respect to 
x. Define the Hamiltonian function by the equation 


fy(x,iiJ,u) = ^gfg (x ,u) + |-?(x,u) 


(5) 


where x and the adjoint variables ip are related by the canonical 
equations 


dip 


x 


( 6 ) 
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if = 


zHfc 

8x 




(7) 


Let u (t) be a candidate control which is admissible and feasible , 

y 

i.e., which satisfies Eqs. (3) and (4). Then for u (t) to also 
be optimal, the following conditions must necessarily be satisfied: 
(a) There exists a number 4 *q <. 0 and an adjoint vector ip (t) , 

such that ip (t) and x (t) are solutions of the canonical equa- 

->* 

tions corresponding to u (t), i.e.. 


X (t) = — ^(x (t), ip(t),u (t);^ ) (8) 

u 

<P (t) = — ^(x (t),^ (t),u (t);ip n ) (9) 

3x U 

satisfying boundary conditions (4) and with a Hamiltonian defined by 
(5). 

(b) The function °H(x* (t),ip (t) ,u ( t ) ; ) is maximized with 

respect to u(t) e ft for u(t) = u*(t), for all t € 

That is, given two admissible controls u (t) and u(t), the 
statements that u (t) is optimal while u(t) £ u*(t) on a finite 
interval of (t^^) imply that 


1ftx*(t),^*(t),5(t);^*) < ^.(x*(t)J*(t),u*(t),^o) (10) 

moreover, ^ = constant 0 and the maximum value, 0, occurs when 
the final time t^ 


is free. 
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8.1.3 Extensions to More General' Optimal Control Problems 

More general cases can easily be imagined, even within the 
class of two-terminal problems. It will be convenient to classify 
them as follows: 

(a) System : (1) autonomous; or (2) non-autonomous ; 

(b) Criterion : (1) integral; or (2) integral plus PVC; 

(c) Terminal Time : (1) fixed; or (2) free; 

(d) Control : (1) variable, u(t); or (2) variable plus fixed 

parameters, u(t) plus 3; and 

(e) Boundary Conditions : (1) fixed or moving target point 

(N-tuple) in R^, x^ or p(t^); or (2) fixed or moving 
manifold (k-fold) in 1^, i.e., g^x.t) = 0, i = 1, 

2 ••• N-k; or (3) free end conditions. 

Necessary conditions, analogous to and derivable from the results 
given above for the basic problems, are summarized here for cases of 
significant interest. For all cases: 

(a) The Hamiltonian is defined as by Eq. (5); 

(b) The state and adjoint variables obey the canonical equa- 
tions (6) and (7); 

(c) The control satisfies the Maximum principle, Eqs . (9) or 
(10). It is assumed here that Eqs. (9) or (10) can be used 
to eliminate explicit appearances of u(t), so that the 
canonical equations depend upon x, jp and t only. This 
may involve, for example, solving the equation 



with 


d % 

du 


3u 


0 


( 11 ) 
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constructing "switching boundaries", etc,; and 
(d) The terminal necessary conditions, applicable to the Hamil- 
tonian and the adjoint variables, are summarized in Table 
5-1 (pp. 306-307) of reference 24. 

It may be observed that the preceding results follow directly from 
the basic theorem. For example, the non- autonomous case is treated 
by simply treating time as an additional state varible, i.e., 

t = x N+1 with x N+1 = 1 (12) 

and then applying the autonomous maximum principle and transversal- 
ity conditions to the resulting N + 1 dimensional system. Also, 
note that the terminal conditions shown in Table 5-1 of reference 24 
follow directly from the partial-derivative interpretation of the 
adjoint variables and Hamiltonian; furthermore, these apply (with 
appropriate changes of sign) to the initial as well as the terminal 
boundary value. 

Necessary conditions for constant design parameters g, • ° * g 

1 L 

may also be derived, by defining L additional state variables 

X N+1 * * ' X N+L Wlth ^+1 = ^+2 * * * " ^+L = 0 (13) 

and again applying the basic theorem. Note that the Hamiltonian, 
and hence also the optimal control law, are unaffected by this step. 
The transversality condition applies to the associated, time-varying 
adjoint variables iJ» N+1 (t) ••• f N+L (t) at both the initial and 

( —V 

final times. If g is unrestricted for example, we have 
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f N+l ft O^ 


- wv ■ 0 


and (14) 

W'f* W'f’ ' 0 

in addition to all of the previously defined conditions. 

8»2 Second-Order Necessary Conditions and Sufficient Conditions for 
a Local Minimum ~ 

Under the present assumptions regarding continuity and convex- 
ity, it can be shown (c.f., refs. 24 and 108) that local sufficiency 
can be established by adding the strengthened Legendre-Clebsch , 
Normality , and Jacobi conditions to the Maximum principle and trans- 
versality conditions . 

8.2.1 Convexity or Legendre-Clebsch Condition 
This states that the second-partial matrix 


§ 


( 1 ) 


(i.e., is negative-definite) for all t in (t Q ,t f ). It simply in- 
sures that is actually maximized, as in ordinary calculus. 

8„2o2 Normality Condition 

A trajectory leading to the point x^ and lying in the inter- 
ior of $tt;t k ,x k ) - c.f., Eq. 8. 1.1(3) - is called normal. In the 
Zermelo problem for instance the two envelope lines for v < u 
(Fig. 3-3) are in fact the boundary of ^(t;^,^). A trajectory 
lying on this boundary clearly does not possess a 2-sided family of 
neighboring extremals; it is impossible to reach x, from any 
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point outside of if (t;t^k x^ in any length of time 
this is manifested by the singularity of the matrix 
3. 2. 2 (7)). Alternatively, normality is verified if 


Mathematically , 
B (c.f., Eqs. 

B is non- 


singular. 

8.2.3 The Jacobi or No-Conjugate-Point Condition 

The Jacobi condition states that an optimal trajectory may not 

have a conj ugate point at which the constant - J contours have dis- 

2 ->-2 

continuous slope, i.e., at which 9 j/9x is unbounded. Regarding 
the initial point as fixed it may be seen that the second partial 
of J, with respect to variations in the terminal point, is given by 
Eqs. 3.2.2 (8) — (9) ; i.e.. 




( 2 ) 


which must be finite. Or with the terminal point fixed, the matrix 


9 2 J 


9x 


= E = 


-B 1 A 


(3) 


must be finite. 

8.2.4 Equivalence of Second-Order Conditions 

The above mentioned normality and Jacobi conditions are often 
seen developed in terms of the backward-sweep approach, c.f., refer- 
ence 24, rather than in terms of the present forward-sweep or 
transition-matrix approach. Therefore, it is of some interest to 
see that the derived matrices E, F, G, and H of Eqs. 4. 2. 2 (9) obey 
exactly the same Riccatti matrix differential equation as do the 
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"Gain Matrices" Q, R, U, and S of the backward-sweep approach of 
reference 24. 

The transition-matrix differential Eq. 4. 2. 2 (3) may be written 
compactly as 


where 


_d_ Ta Bj ra el I A B| 

dt !_c dJ Ly §J L c dj 


— f a 9g j 9g 

a - — , B - — , y = and 6 = — 

9x Sip 3x 3$ 


(4) 


(5) 


By formally differentiating the defining Eq. 4. 2. 2 (9) for the de- 
rived matrices and using the convenient identity that 


(B 1 ) = - B 1 BB ~ 1 


( 6 ) 


i t is readily computed that 


(7) 


E = FgG 
F = FBH - Fa 
G = HgG - <5G 
H = HBH - 6H - Ha - y 

These are Riccatti differential equations, identical in form to Eqs 
5. 3. 3. 4, 5. 3. 3. 5, 5. 3. 3. 8 and 5. 3*3. 9 of reference 24. 

The latter, in the present a, 8, y» and S notation, may be 
written: 

Q = uBR 
u = uBS - ua 
R = SBR - 5R 
S = SBS - <5S - Sa - y 


( 8 ) 
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Equations (7) and (8) do not possess finite initial conditions and 
hence, cannot be integrated forward along with the canonical equa- 
tions. Applicable terminal conditions, however, are given in Chap- 
ters 5 and 6 of reference 24. For example, if the terminal bound- 
ary condition required is to attain a single fixed point 

x(t f ) = x f 

then 

Q(t f ) = 0 
U(t ) = I 

(9) 

R(t f ) = I 
S(t f ) = 0 

Because of the symmetry of Eqs . (8) and 4. 2. 2 (4) and the form of 
Eq. (7) it is clear that 

U = R T and G = F T (10) 

Thus, the present derived matrices would obey the same Ricatti dif- 
ferential equations as those arising in Hamilton-Jacobi theory, but 
the boundary conditions are different. 

The precise relation between the present E, F, G, and H 

matrices, and Q, u, R, and S of reference 24, may however be dis- 

2 

played by considering the second variation 6 J. In terms of the 
present analysis we have that 


6 J = [<5x Q ,6x f ] 


E ] F 


6x n 

Gji 


1 

1 X 

I 40 ! 


( 11 ) 


while according to page 183 of reference 24 it is 



( 12 ) 


5 



[<Sx 0 ,Sx f ] 


-IT 1 T -1 



S-RQ R ! RQ 


6x q 

-1 I -1 



Q R | -Q 


ox^ 


By comparing (11) and (12) the following relationships evidently 
hold: 


E - -B 1 A 


= (S - RQ V)| 
c 0’ 1 


|t.t r 


F = B 


-1 




T -1 

R Q 


t,t 


0 


G = ( C - DB 1 A) = Q 1 R 

V* 


t»t f 


H = DB 


.-1 




t,t r 


From this , it is clear that the "normality" and Jacobi conditions as 
developed in reference 24 using the "backward-sweep" approach,' are 
exactly equivalent to those given here* 

8.3 Convergence of the Transition Matrix Algorithm 

A schematic of the algorithm is shown in Figure 8-1 below. 

Its convergence is demonstrated as follows. Let £ = x(t^) - x^; 
this depends on |(t Q ) only. That is, £ = £(iL) . To find the 
roots of £( 4 > q ) = 0, just expand 

toL) = ?oL) + 

i+1 i Usl 

or, dropping subscript 0 we have 


3|(l 0 ) 


->■ 

3'ifr 


(i) 
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Perturbation 

matrices 


FIGURE 8-1. - SCHEMATIC OF TRANSFER MATRIX ALGORITHM. 
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* 1+1 


*1 " 


- ■+■ 
85 


l4. 


-i 


-+ 

5, 


( 2 ) 


Thus 


■f ->* 

*1+1 " * = *1 " 


Hal 1 1 


■>* 
? i - * 


r+T 

*i 


95 


9^ 


-+ 

5, 


h-* 

- u - 


■+ 

95 


L94d 


-+ ->& \ ->- -i. -V3&T 

5(^ )j = <K* ± ) - <P(t ) 


(3) 


Next, expand <J> in Taylor's formula around For the j th com- 

ponent we have 


4- -> ->* ->* T 

= 4>j Cip ) + ('P i ~ <P ) 


thus 


— — I 

9 j> . (<P ) 


9^ 


1 -*-* T 

+ y 0P ± ~1> ) 


2 

9 *.( 0 ) 

*L 


a| 2 


•>* 

01^ -4> ) 


( 4 ) 


-> ->* ->- ->- A T 9^4 1 _i_ rp 

(^ i+ l ~*P ) . = (<f> i , =('P i ~IP ) r — J+Y 0p ± ~lp ) T 

1 J 3 ip 


2~> ->■ 

a 4>. (0) 

— J 


L 


->• ->A 

(* ± ) 


(5) 


->* 


where Q - + (1 - Q)ip and 0 < 0 < 1, But, from (3), 


* 

9 (j> _ 

dip 


[I] - 

•+ 

95 

-1 

(— 

95 

->* 

+ 

H 

-+* 


_3lJ> _ 


dip 


_dlp_ 


-2 2-y 

hkm -o 

9i|j 


( 6 ) 


Thus, Eq. (5) leads to the result that for i > some critical value 
we have quadratic convergence , i . e . , 


-f* i 2 


* i+ l " ** I -2- P I ~ 'P 


( 7 ) 


where 
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P 


max 

j 


max 

e*[o,i] 


r 

2 


J max +T 

3 4> ± (0C0)) 

v > 

|y|#> y 

w 

+2 

dip 

y 


( 8 ) 


Directions of descent . - It should be understood that the quad- 
ratic convergence displayed above is an ultimate rate, attainable 
only in some "sufficiently Close" neighborhood of the actual root, 

• Thus, initial rates of convergence are of equal concern, be- 
cause it is evidently possible for a poor choice of to cause 

the Newton algorithm to diverge. Fortunately, the step direction s 
defined by Eq. (2) is also a direction of descent for the quadratic 
scalar terminal error function. 


I ll(?)l 2 

That is , since 


we have that 




2 = f - 1 it 
^ -»• 
dip 


(9) 


i 1 2 ■+ -*T 

Vj|5| •> - - E 


H 


dip 


H 


dip 


-1 


->• 

Z 


- Ill 2 < 0 


Therefore, by the theorem in Section 3.2 , convergence, and hence 
by implication the quadratic convergence domain can eventually be at- 
tained as long as successive directions of descent are well defined. 
These directions are well defined if the matrices dt(ip n )/dip are 
non-singular. The singular case need not stop the descent if the 
gradient vector shown in Eq. (9) is non-zero since then an alterna- 
tive direction of descent is defined. Even at a spurious (non-zero) 
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local minimum, it may be possible to replace the error function \%\^ 
by a more general form, such as £ Q£, where Q is some non-negative 
matrix, and thereby recover a legitimate direction of descent. 

Thus, it may be seem that under fairly general circumstances, 
the convergence rate of the present algorithm need be no worse than 
linear, initially, and will ultimately improve to quadratic. More- 
over, it is generally possible to have quadratic convergence from the 
very beginning if the interval of integration At = is 

small enough. This may be seen as follows. For the first step, 

Eq. (7) may be written as 



2 


where the spectral norm. 


max 



max 


y|#> 


y L 


9 ^( 6 ( 0 )) ^ 
y y 


-+2 

dip 


(8a) 


is the same as the spectral radius of 


3 2 <f> . 

l 

->• . 
3x3tj; 


since this matrix is 


symmetric. It is well known in algebra that the spectral radius of 
a matrix is bounded by row and column sums , viz 


( 10 ) 



n 

2 

n 

\ \. 

2 

<. min<; 

max \ 

3 <J> i 

max \ 

3 <j>. 

1 J Lj 

3^ 3^ 

5 k LJ 

9^9^ k 


k=l 


1=1 


Referring to the definition of <J> it may be seen that 
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3* 2 



c 2 

dip 


3ii 

->3 

9<r 


^ v_2 /7 1 

Va|' 2 


9i|) 


3£> 73~g 


4 1 

££ 

9^ 


(ID 


Because of the definition of £, it follows that — = — and simi- 

dip dip 

larly for the higher derivatives. Equations (3a) and (4a) then show 
that, for short times, 

9x(t) _ 9? 


9^(t 0 ) 


i • At = <y(At) 


dip 


( 12 ) 


Also, 


& 

dip 2 


> 2 t 


9f 


dlpdx dip 


At = ©'(At) 


and in a similar fashion 


„3-> , 

■“§ = Q^(At 3 ) 

dip 

Thus, each term in Eq. (11) is of order O'(At) and, hence vanishes 
as At 0. Therefore, p, which is bounded above by sums of such 
terms, is also of order ©((At). This means that we can make p as 
small as we please by choosing a small enough At; the criterion for 
quadratic convergence can be satisfied in this manner for an arbi- 
trary ip ^ as long as d^/dip is not singular. 



9. APPENDIX B 

DISCUSSION OF LEVEL II - SUCCESSIVE IMPROVEMENT STEPS IN STATE SPACE 
At this stage, the general problem given in Section 2.1 has been 
divided into a series of "short" sub-problems by imposing mesh 
points, and solutions to these short sub-ptoblems have presumably 
been attained. It has been shown that J depends on mesh point co- 
ordinates only, the gradient of J in X has been derived, and 

pp 

the "Necessary Conditions" for unconstrained mesh points and for 
simple constraints have been developed. 

9 . 1 Kuhn-Tucker Necessary Conditions 

The celebrated Kuhn-Tucker conditions would apply in the pres- 
ence of more general constraints on X . See pages 17-34 of refer- 

PP 

ence 45 for a complete discussion. Let the equality-constraint man- 
ifolds be described as follows : 

! 

% = ^|h jk (x) = 0, j = 1, i^} (1) 

where 

k = 0, 1 • • • K + i (2) 

Similarly, the inequality constraint may be described as 

r = {x|g i (x) 0, i = 1 ••• m} (3) 
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Briefly, the first order Kuhn-Tucker condition states that J is 
minimized when no feasible perturbation of X will bring about a 
lower value. 

Geometrically, this means that V„J must lie within the cone 

xL 

formed by (a) the gradients to the equality-constraint surfaces, 
V X h jk®’ and ^ the § radients to the binding inequalities, 

^X®b ’ ^ w ^ ere beB, the set of values of all binding constraints) 
Algebraically, for V^J to lie within the above described cone 
it must be expressible as a linear combination of the gradients to 
the equality and binding inequality surfaces, i.e., 

K+l n k 

V 1 " L, X b V X 8 b (X) " / / ) ,/ u jk V X h jk (X) 
beB jc=0 j=l 

Here and are "generalized" Lagrange multipliers, respec- 

tively associated with the inequality and equality relations. The 

I 

inequality multipliers satisfy the further condition that 

^ 2. 0 if g (X) = 0 (constraint binding) (5) 

and 

^ = 0 if g^(X) > 0 (not binding) (6) 

Conditions (1) to (6) above, together with a certain regularity as- 
sumption (the so-called First Order Constraint Qualification, c.f., 
p. 19 of ref. 45) comprise the Kuhn-Tucker Necessity Theorem. 

9.2 Descent Via First-Order Gradient Steps 

The "gradient" techniques described (e.g.) in reference 81 have 
the desirable traits of being conceptually simple, straight forward, 
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and easy to apply. Unfortunately, initial numerical examples (c.f.. 
Section 4.1) displayed unexpectedly poor initial rates of conver- 
gence, especially when a large number of mesh points was used. In- 
vestigation showed that this is due to the extreme slopes or curva- 
tures, which tend to develop near a fixed point, having a dispropor- 
tionate effect on the value of the criterion J. This, in turn 
sharply limits the step size that can be taken. 

These considerations suggest that descent directions which van- 
ish at the fixed points and which are in some sense as "smooth" as 
possible, may yield better computational results. Such a direction 
would allow points far from a fixed terminal, to be moved through a 
relatively large displacement, while those closer to a fixed point 
are limited to small displacements. In the next section, we derive 
a direction s, having this quality, as a solution of a variational 
problem. 

9.2.1 Auxiliary Variational Problem 

Minimize 


n t f 


fc o 


s ( n > 




dt where s 


(n) 


A|) 


dt 


n 


( 1 ) 


subject to the boundary conditions 


' * s(t jL ) = 0, i = 0,1,2, •*' K 
and the isoperimetric constraint 



( 2 ) 


S ‘ Aip = constant < 0 


(3) 
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Equation (1) guarantees that s will be maximally smooth to order 
n; (2) results in the automatic satisfaction of the fixed point 
boundary values, and (3) guarantees s to be a direction of descent 
and provides a common basis for comparing different directions. (Al- 
though the problem is formulated for continuous time, with Arp(t) 
as defined as in Section 3.2, the formulation and results also apply 
to the discrete case, i.e., by defining Aip for Eq. (3) to to be an 
impulse function “ 5(t - t^)) . Introducing the multiplier A and 
writing 

F(s,s (n) ) = |s (n) | + 2As • A^ (4) 

the Euler-Lagrange equations for this problem become 


n 



(-l) j 




= 0 


(5) 


and, in case n > 1, the "natural boundary conditions" (refs. 53 and 
173) 






= 0 


( 6 ) 


apply for j = 0, 1, •*•, n - 2. In terms of the problem at hand, 
conditions (5) and (6) become 
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(-1) s = - AAip 


s(t Q ) = s(t f ) = 0 


i M (t 0 ) v - o 


(7) 


Kn+1) 


(t„) - J (n+1) ( tj - 


J (2 "- 2) (t 0 ) .?< 2 "- 2 >(t f ) -0 


Thus, we finally obtain s as a polynominal in t, augmented by a 
2n-fold quadrature. 

It is of interest to verify that the step s thus defined is 
indeed a direction of descent. Consider the inner product 


= A 


s • Aip dt 


'0 

)t. 


= A 


-*■ ->-( 20 ) , i n 

s • s' (-1) dt 


= a Ht • J^- 1 ) 


tf _ -(l)-(2n-2) ,tf 
t 0 "0 


+ J (2) ! (2n - 3 > 


+ (-D 


n 


:(n) 


s (n) dt> 


( 8 ) 
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Here again the boundary terms vanish, because of (21) and the inte- 
grand is positive definite wherever Ai|> ^ 0. Hence, by appropriate 

choice of the sign of A we can always achieve that p < 0 (if 

A 

A ’I' ^0), i.e., s is a direction of descent. Therefore, the theorem 
in Section 3.2.3 implies the eventual convergence of the following. 

9.2.2 Linear Step Algorithm 

(a) Xg(t) = feasible, but otherwise arbitrary initial guess 

(b) Ai^(t) = as defined in Section 3.2.2 

(c) Sj(t) = solution of Eqs . (5) to (7) above with |a| =1 

(d) A = chosen to minimize J(x (t) + As.(t)) 

J J J 

(e) x j+1 (t) = x^(t) + *jSj(t) 

(f) Pj = as defined by Eq. (8) above 

(g) Iterate (b) - (f) until the estimated decrement reaches a 

satisfactorily small magnitude, i. e., 0 |p.| < e 

9.2.3 Ultimate Convergence Rates 

Under the conditions applicable in this section, the theorem in 


Section 3.2.3 implies the ultimate convergence of the preceding al- 
gorithm. It remains to examine its rate of ultimate convergence. To 
do this, note that the recurrence relation s£ + ^= - | A | A^ k may be 
written as 


£< 2n) 

S k+1 





+ As, ) + 
k 



(j+1) 


+ 


3h 

at 



+ As k ) ,£(x k 


Xs k’ x k 


+J 

As, 

k 



( 1 ) 
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That is 
depends on 


s (2n) 

’ S k+1 

s k’ s k’ 


may be written as a composite function 

. , ->(J+1) , 

• s, alone. That is 

k 


-y 

q 


which 


s 


(2n) 

k+1 


/ * •• 

q \VW 


^(J) +(J+1)\ 

S k ,s k y 


( 2 ) 


This shows that the present algorithm is formally equivalent to 

Picard iteration. The convergence rate may be exhibited by subtrac- 

til 

ting the k such equation from the k + l 1 " : 


g" ( 2n ) _ +(2n) _ 
k+1 s k 






(3) 


If it is now assumed that the composite function q is continuously 
differentiable with respect to its arguments ~s, s, •••• s^ J+1 \ the 
mean value theorem for derivatives may be used to show that 


+(2n) _ +(2n) = 9<L . _ + ) + + n 

k +1 S k + l k k-l' + (s k s k-r 

dS dS 


+ 


ii. 


(j+i) 


8s 


r(J+D 


HJ+D) 

s k-i y 


(4) 


Now define 


-> -v -y 

V k+1 = S k+1 " S k 


(5) 


then 


->-(2n) 9q -> 8q 4- 

v k+l - • v k + Tt ' v k + 

OS 8s 


+ 


*3 V ( J +D 


8s 


(J+l) k 


( 6 ) 


Now 2n integrations lead to 
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v. 


k+1 


r 

/ • • • i 

I 5i * d t 2n + / 

r / 

/ oJ 

£ k J 



■5a t dt 2n + 

d k 


Next, applying the mean value thus for integrals yields 


(7) 


v = Isl 

k+1 + 

dS 


-* , . 2n , 9q 

v. dt H ?■ 

k 

0“ 8s O' 


v, dt 2n + 
k 


= ia 

9s 


2n 3f 


?k dt "’ + d 

0 v ds ~ O' 


j*.2n-l 

v, dt 
k 


+ -4! 


9s (J+1) 


£ dt 211 - 1 " 1 


(8) 


Finally, by taking magnitudes, using the triangle inequality and 
properties of the scalar product, and applying the maximum operation 
yields the result that 


max|v k+1 |<. max 


b 2n ^ 

2 ^t max | v k | + max 


iSL 

9 s 


9£ 

9s 


b 2n “l y 

"(2n-l) I max ^ k l + 


• + max 


Ja. 


9s J+1 


b 2 n- J -i 
(2n-J-l) ! 


(9) 


That is. 


J+l 


max l^k+ll < max I^J 

t ‘S'X t e'ZC 


\ 

J 


max 


.laL 


9s 


U) 


£=0 


b 2n £ 
(2n-£) ! 


( 10 ) 
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Therefore, assuming that 2n >_ J + 1 and that q is sufficiently 
differentiable, it has been shown that the present algorithm will 
converge at least geometrically if the interval b = t^ - t^ is 
chosen small enough. The result may be extended to longer intervals 
by partitioning into I sub-intervals and then extending the maximum 
operation to include the stipulation that i <. 1 I. 

9.3 The Hessian Matrix and Descent Via Second-Order Steps 

After having selected a set of mesh points x^ and performing 
the point-point transfers as in Section 3.1, the following informa- 
tion is available at the beginning and end of each arc: 

(a) The values of |i(t) and A^(t) 

(b) The matrices A = B = C = and 

3x(t k ) 3$(t k ) 9x(t k ) 

D „ J?(t) 

9^(t k ) 


(c) And of course the values of x(t) itself 
This same data that was generated in the process of computing the in- 
dividual two-point transfers can also be used to define a second- 
order improvement step, y(t) , with the following properties: 

(a) The direction 


Y = [y. 


y K ] ' 


is a direction of descent in R 


KN 


(b) The sequence x^(t) + YjCt) converges quadratically to 

-V 

. (t) for sufficiently large j and if 


x 


opt 


max 

t 


V e) - *opt (t) 


is suitably small. 
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9.3.1 The Hessian 

As in Section 3,1, the perturbation equations may be written 


y 


"A 

B 


y 

1 

k+l 

_C 

D 

k 

-V 

-l 


( 1 ) 


or 


where 


\ ' 


E F 



\+l 


G H 

k 

ik+i 


( 2 ) 


\ 


\ - \ 1 


G k = [C - DB A] k and = . D k „ k 


-1 


(3) 


Applying this at the beginning and end of each subarc we see that, 
for the x Q - x^ transfer, 


*0 ■ Vo + Vi 


? 1 * G 0?0 + Vl 

Similarly, for the x^ - transfer, 

?' = E 1 ? 1 + Fj_5 2 


*2 = G in * v 2 


(4a) 


(4b) 


after which, 


^2 = E 2 y 2 + F 2 y 3 

+ + _ T>- 

*3 = G 2 y 2 + H 2 y 3 


(4c) 


and, so forth, until finally, for the K to K + 1 transfer we have 



- 157 - 


E K y K + F K y K+l 
^K+l = G k y K + H K y K+l 


(4d) 


To first order in small quantities, the perturbed gradient ¥ in 
R^j may be expressed as 

' Atj^ + A^ 


¥ + AT = 


+ a K 


(5) 


Therefore, y must be chosen so that A^ = - AiL, 1 <. k < K. By 
subtracting Eqs. (4a) from (4b), (4c) from (4d) , etc., it is seen 
that 

- 4 *i ■ - *1 - = - G 0 F 0 + (E 1 - V*1 + F 1?2 


( 6 ) 




- 4 *K - 4 *K - +K - *K - - G K-1?K-1 + (E K - VA + V K+ 1 


• *T 

Since y Q = y K+1 = 0 by definition, this set of equations may be 


written in matrix form as 



K 


Aip. 

A^>, 


A\p, 


K 


( 7 ) 
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or, more compactly 


KY = - 'V 


( 8 ) 


where K is the KN x KN, block-tridiagonal Hessian matrix indica- 
ted above and Y and ¥ are as previously defined. Thus, if K 
is non-singular the solution 

Y = - K -1 ¥ (9) 

may in principle be computed immediately since the A^’s and the 
elements of K are defined in terms of the preceding trajectory. 

9.3.2 Quadratic Step Algorithm 

(a) The initial mesh, x^ q, is arbitrary, but feasible 

(b) The mesh x^ ^ is defined at the beginning of iteration 
number i. Compute the optimal sub-transfers per Level I. 

(c) Using the AiJj's and sub-matrices thus defined, form 

K, and then solve for the mesh-point perturbation vector 
Y 

(d) Convergence is attained if |Y[ is sufficiently small. If 
so, stop; otherwise continue 

(e) Define the i + l*"* 1 mesh as 

J k,i+i ' *k,i ’ *k,i 1 * k * K 

and also estimate initial Y's for the next set of 2-point 

sub-problems, i.e., i£(t k ) = |(t k ) ■ + £“| , 0 <_ k <. K + 1 

i+1 i K i 

by evaluating Eq. (2). Repeat (b) - (e) as necessary. 

9.3.3 Convergence 

Ultimate convergence may be most easily studied in terms of an 
auxiliary function, which is here taken as | T [ . The gradient of 
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this function is 

\M 2 = * T [f] = * t k (d 

But the search direction is Y = - K - '*''?, hence 

V X I ^ I 2 * Y = ~ ^ T K • K _1 ¥ <. 0 (2) 

9 

so that Y is indeed a direction of descent for | ’{' | . Thus, ultim- 

ate convergence (to |w| » 0) is implied by Section 3.2.3 provided 

that K and its component parts are well defined and non-singular. 

2 

Indeed, it can be shown that f | -> 0 quadratically for n 
large enough. Consider the recursion formula: 


n+1 


_ «, (31] 1 

n 

' 'n 


Y N ■ ’n - K A - *<V 


Expand $(¥ ) in Maclaurin series - 
n 


M + 1 »T I 

a 2 < j> C © ) 


2 

L J ¥=0 | 

_ a’f J 


n 


But 


9£ 

dV 


= 0 


^=0 


Hence 


V . , = V 
n+1 n 


9 2 $(6) 

dV 2 


n 


or, again we have quadratic convergence 


n+1 


<. P I'V, 


n 


(3) 


(4) 


(5) 


( 6 ) 


(7) 


where 
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1 3 V (0) 

l 


ar 


( 8 ) 


The convergence of | T | to zero could imply a maximum as well as a 
minimum, value of J. For a minimum to occur, the matrix 


K = 


3T 


3X 

must be positive definite, and then the estimated decrement in J is 


AJ = Y 


¥ = - Y KY 


( 10 ) 


as long as | Y j ^ 0. In other words, convergence is not only quadra- 
tic, but monotonic. 

9.4 Descent Via Fletcher-Powell Algorithm 

This may be accomplished, for example, with a modified Fletcher- 

Powell routine. That is, successive directions of descent s are 

n 

defined in 31 by the formula 

pp 


Vi ■ - o> 

where L n is an appropriately dimensioned positive matrix defined 
by the recursion 


L , = L + 


-*■ -+T -+■ ->• T 

a • a (L <f> ) (L cf. ) 
n n n T n n Y n 


n+1 n ->T-> 1 T, t A , 

o ( p (p (L d) ) 

n Y n Y n v n y n ; 


( 2 ) 


where the change in position is 


-+->■-> 
a = x - x , 
n n n-1 


and the change in the gradient is 


“>■ 

<p = Arp - Arp , 

n r n r n-l 
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The successive step-lengths, ct n> are predicted from the second- 
partial data instead of being determined by unidimensional search- 
ing. Expanding J n+ 2 to sec °nd order in 6x yields 


t>t 1 ->T -*■ 

J * J + Ai/j <5x + — 6x K 6x 
n+1 n r n 2 n n n 


(3) 


Replacing 6 x r by a n s n > equating 8J/3a n to zero and solving for 
a n leads to a predicted step length 

A \p s 

n n 

“n " " ->T -► 

s K s 
n n n 


(4) 


which (to second order) will minimize J along the s -direction. 

n+1 n 

To summarize, the preceding descent algorithm in consists of 

the following steps : 

(a) Lq = arbitrary positive definite matrix, e.g., Lq = I 

(b) J = - L a| , per (1) 

(c) a n per (4) , or by numerical search if (4) does not lead to 
a reduction of the criterion J. 

(d) x = x , + a ,s , 

n n-1 n-1 n-1 

(e) L n+1 per (2) 

Steps (b) - (e) are repeated until an appropriate measure of conver- 
gence has been satisfied. 

It can be shown that, for this algorithm, 

(a) Lim L = K 

n n 

n->°° 

(b) Lim a =1; and 

n 

n-*» 

(c) Convergence is quadratic. 1 

See references 48, 49, and 68 for further discussion of this algorithm. 



10. APPENDIX C 


AN ANALYTICAL APPROXIMATION TO THE FUNDAMENTAL MATRIX FOR 

SPACE TRAJECTORIES "" ' ' ' : ' ' 

For t;he purpose of this section, let the vectors v^, , v^. 


and denote the unit solutions of Eq. 4.4(10); i.e., they are 

the column vectors that comprise S>, and their individual component 
equations for column i (where 1 <. i < 4) are 


4 - -> 

V li " V 2i 


v 2± = - to (t)^, +' 


li 


4i 


3i = “ (t)v 4i 


V 4i V 3i 


> 


with 


In second order form: 


v ji ( V 6 ji 


v 4i = " “ (t)v 4i 


( 1 ) 


( 2 ) 


li 


= - 0) (t)v Xi + gv 4i 


(3) 


If to were constant, the above equations would yield harmonic solu— 
ions. These, however, are quite inaccurate except for trajectories 
involving only a small variation of radius (and hence, to). 
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On the other hand, Eqs. (2) and (3) can be solved analytically 
for a special, but more general, form of tp(t) which can be made to 
exactly match the "real" function at the beginning and end (or any 
other two points) of the subarc. 

That is , let 


JL _ v 
2 

(0 


a constant 


(4) 


hence 


'k+1 


= Kt 


k+1 


or, assuming for convenience that t^ = 0, to(P) 


“Q, and w(t k+1 ) 


we find that 


K = 


a) 0 t k+l a) f t k+l 


or 




where t ' == t - t is the elapsed time since the beginning of the 
sub arc . 

A change of variable . - Now consider the transformation defined 
• • ... 
by t' ■> z, where z » w. Hence, 


d( ) _ 
dt 


0) 


d( ) 
dz 


( 6 ) 


i 2 ( ) 2 d 2 ( ) . dm d( ) 

— -7T~ - 0) — 7T- + tO — T 

dt 2 dz 2 dz d? 


and 
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Thus, the substitution t -* z reduces (2) to the form 


v" + ui'/o) v! . + v. . = 0 
4i 4i 4i 


( 7 ) 


I • r\ 

But, since to /<o = to/ur = K where K is a constant (c.f., Eqs. (4) 
and (5)) Eq. (7) reduces to harmonic form 

( 8 ) 


v 4i + Kv 4i + v 4i * 0 


The homogeneous solution . - The solution of (8) is well known 
to be either 


Kz 


if K <. 2, or 


= e [A sin x + cos x] 


_ M 

2 r 

= e sinh x + cosh x] 


if 


K| >2 


(9a) 


(9b) 


where 


x = 


K 2 

1 “ T 


1/2 


z = az 


( 10 ) 


The oscillatory solution (9a) will be seen to be of primary interest. 
Since K was defined as 

1 1 


K = 


(0_ t ' (0_t ' 

Of f f 


(where t^ = t^ - t^) it is clear that |k| will, in general, be 
small if t^ is large. That is, if 
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(U) 


we need never consider (9b). For trips from the Earth inwards, 
o)q = 1 and > 1, so that (28) is always satisfied easily. That 
is, we need t^ > of the Earth's Schuler period (about 29 days) 
if we are going tq 0 A.U.'s (w^ «) . Fof outward transfers, fhe 

final angular velpcity can become very small, and we must con-?- 

sider that t f increases in proportion tq the period of the final 
orbit. This does not seem an unreasonable assumption -r howevsr, 
condition (11) should be verified numerically when "fast" outer- 
planet trips are being considered. But, for the present we will 
only consider solution (9a). 


The transformed argument , - In order to work with Eq. (9), in 
either of its forms, we need to know z as a funqtfon of f and 
a) as a function of z (as well as t) . first, let us integrate 


Eq. (6), i.e-» 


i.e. , 



(12 a) 


or 



(12b) 
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The Inhomogeneous solution . - By applying the same transforma- 
tion to Eq. (3) we obtain 

v ii + Kv ii + v h * -2 \i 

0) 

where is the homogeneous solution given in (9) above. First, 

the complimentary solution of (14) is of the same form as (9) namely 


Vu(z) = exp [C i sin x + D. cos x] 

complimentary 

Using the method of variation of parameters we try a particular 


(15a) 


solution of the form 

v 1;L (z) = exp [U(z) sin t + V (z) cos t] (15b) 

After some tedious but straight-forward calculations it is found that 
U(z) and V(z) must satisfy the auxiliary differential equations 


V| sin x + cos x = 0 


v ; 


a cos t - -j sin t 


rj - Uj ji 


a sin t + y cos x 


1 A i 

r J" 7 


B, 


sin x + — j cos x 

CO (0 


(16) 


which leads to the result that 


V*> --T ^i( e ' 2KZ 

a “o V 


-K sin 2x -! a cos 2x 


2 2 

4 or + a Z ) 


+ B, 


L-2KZ 

-2K cos x +2a sin r 

V ae' 2Kz 

\ I 

2 2 
L 4 (K + a-) 

/ 4(K 2 a 2 )j 


( 17 ) 
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and 


U(z) = — ^ K A 


am 


0 


/ -2Kz 


-2Kz 


2 2 
4(K + a ) 


+ B. e 

l 


-Kz 


[- K cos 2Y + a sin 2 t] 


K sin 2x - a cos 2x 


2 2 
4(K Z - a 4 ) 


(18) 


Then, putting (17) and (18) back into (15b), adding the compliment- 
ary solution (15a) and simplifying we finally obtain 


5Kz 

2 


V li . 2 r 2 ^ 2. 

4Wq[K + a ] 


h 1 


a 


Aj, sin t + — cos x 


t J + 


cos t - — sm 
K 


Kz 

2 


sin t + D i cos t; 


(19) 


This together with the homogeneous solution (9) represents the gen- 
eral solution for one column of 0. 

The typical column may be written as 


with the constants A. 


so that 


v i = 


V li 


mv 


li 


" v 4i w 


L V 4i 


( 20 ) 


D ± chosen so that $(t Q ;t 0 ) = I; i.e.. 


VV ' V 


( 21 ) 
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Define auxiliary functions f^ f 2 , . . . f ; g ;..g 4 and 


follows ■ 


•h, as 
4 


5Kz 

2 


f l , 2 r „2 . 2, 


4« n [K + a ] 


_ -5Ke 


5Kz 

2 


f = 

1 ~ 2 r „2 


f„ = 


8^tr + a 2 ] 


sin t + — cos t 


f 2 


f. = 


a 

a cos x - — sin x 
K 


cos x - — sin j 
K 


f 3 


> 


-a sin x - — ■ cos x 
K 


f , = e 


Kz 

2 


f i 


-K 


Kz 

2 


f- = sin x 


f 5 


= a cos x 


f, = cos x 


*6 


= -a sin t 


(22) 


Si 


So = 


So 


S/, 


f l f 2 

f l f 3 

f 4 f 5 

f 4 f 6 


( 23 ) 
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h i - 8i - + f{f 2 


h 2 g 2 f l f 3 + f l f 3 


> 


Using these functions may be written as 


v i = 


*li 


Sl A i 

g 2 B i 

g 3 C i 

g 4 D i" 

"*li 


(oh.. A. 

1 l 

0)h o B . 
3 l 

mh„C. 
3 l 

wh.D. 
4 l 

-“+41 


-mh-A. 
3 l 

-rnh . B , 
4 i 

0 

i 

0 

j 

•ri 

<3* 

^ ! 


g 3 A i 

g 4 B i 

0 

0 


which may be solved to yield 


A 


g i 

g 2 

g 3 

•«1 

-1 

1 

< 

I — 1 

h- | 

B 



wh 2 

wh 3 

toh , 
4 


V 2i 

C 


-cohg 

-uh, 

4 

0 

0 


V 3i 

D 

i 

_ g 3 

g 4 

0 

0 


_ V 4i 


provided that the inverse indicated in Eq. (26) exists, 


( 24 ) 


(25) 


(26) 


That the inverse does exist is readily shown :by computing the de- 
terminant ; i . e . , 


A - - 


g 3 

g 4 

h 3 

h 4 


-ooh„ -toh. 
3 4 


f 4 (f 5 f 6 - f 5 f 6> 


2 -2Kz 2 , n 

= - a e = -a when z = 0 


Since a^ = 1 - K^/4 we may rely on the same assumption ( J K | < 2) by 
which we selected the oscillatory solution (20a) for further analysis 
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The general algebraic form of (25) evaluated at the Initial 
point z = 0, may be written as 

1 


li 


6 


2i 

= 

6 


3i 



4Ku) 2 [K 2 + a 2 ] 1 4 Ku 2 [K 2 + a 2 ] 1 I 

------ l r- 

I _ 5 a_ I 

~3au) | ^ ~ 2 K ~ K | ' -Kio 

22 2 22 2 aco i9 

8a)Q[r + sl] | 4u)q[IT + a Z ] \ 1 1 

1 J__ 


-aoj 


Kto 

2 


-t 


4i 


The solution of (27) is 


h 


(27) 


B i = 6 4i 

A i = 2^ 6 4i “ ^ 6 3i 

1, 'I . 

2 °4i ~ Kui °3i 

1 = 11 " T" 2 'rv 2 * 2l 

4(*)q[K + a ] 


(28a) 


and 


C i 2aoj 6 2i + 2a 6 li + , 2 r 2 


4au)g[K + a ] 


26 


3i 


- (I k - 4 + 1) «« 


It is also of interest to compute the values of A. ••••D. that qre 
appropriate for an impulsive- thrust solution (i.e., 6=0 in Eq. 
(3)). Following the procedure shown above, it is readily seen that, 
when 8 =0, we have 

A i = 6 4i “ it 6 3i 

B i =6 4i 

r 

P _ 2i , K 6. . 
i 2aw 2a 1 


D i = 6 li 


(28b) 
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Analytical form of the fundamental matrix . - We have now seen 
that each of the 4 column vectors v . (i = l.,.,4) composing 
3 >(z;Zq) may be written as 



V 


g l 

g 2 

g 3 

S 4 


V 


V 2i 


oih^ 

wh 2 

tuh^ 

o)h. 

4 


B i 


V 3i 


-uh^ 

-a)h. 

4 

0 

0 


h 


_ V 4i_ 


g 3 

g 4 

0 

0 


D. , 

X 

constants 

A i* 

. . .D . 
1 

are 

evaluated by ap] 


(29a) 


(i=i. . .4) 


(22) through (28) successively for i = l.,,4. The auxiliary func- 
tions gj and hj (j = 1...4) are determined by evaluating Eqs . 

(23) and (24) at the time of interest. That is. 


$(z;z Q ) = 


g^' 2 ) g 2 (z) g 3 (z) g 4 (z) 


A 1 A 2 *3 A 4 

^h ( z ) uh 2 (z) ooh 3 (z) cuh 4 (z) 


8 1 B 2 B 3 B 4 

-ti)h^(z) -a)h 4 (z) 0 0 


C 1 C 2 C 3 C 4 

g 3 (?) §4(2) 0 P 


D 1 D 2 D 3 B 4 


Or, renaming the above sets of functions and constants as 


(29b) 


G ■ [ s lk i 


and 


c ' 'V 


we may write 
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$(z;zq) = GC 

so that 

V 2;z o> = B vvi 

Summary of subarc computations . - (a) Compute the auxiliary 
function m(t) and K (Eq. (5)); z and u>(z)(Eq. (12)); a and t 
(E q. (10)). (b) Evaluate Eqs. (22) through (29) at z = 0 to find 

the 16 constants A^...D^ (i = 1,...4) and the 16 components of the 
fundamental matrix $. (c) Partition $ onto the A, B, C, and D 

blocks as in Chapter 3 and, using the given value of 3 together 
with the elements of these blocks, solve Eqs. 4.4.1(15) and 
4.4.1(13). (d) Equation 4.4.1(11) may be evaluated at intermediate 

times, t^ <. t <. t^ + 2 to obtain a picture of the state and adjoint 
trajectory, (e) Numerically compute 



Analytical computation of J . - As an alternative to numerical 
integration, it is possible to directly compute that (e.g.) 


where 


y l (z) = A i§ 3 (z) + B i® 4^ z ^ 


A 1 A l x 0 + A 2 v 0 + A 3 X 1,0 + A 4 y l,0 


(30) 


B 1 B 1 X 0 + B 2 V 0 + B 3 A 1,0 + B 4 y l,0 


li 


( 31 ) 
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There are similar expressions for y 2 a P<i P 3 • Note, however that 
although the same functions § 3 ( 2 ) and g^(z) are used for the y's, 
the constants (31) are in general different. Hence, let us rewrite 


(30) and (31) as 


Pj (z) = Ajg 3 (z) + B_.(z)g 4 (z) 


where 


and 


Thus , 


A_. = A^Xj (0) + A 2 Vj(0) + A^lj (0) + A 4 y_j(0) 


B. = B..X . (0) + B 0 v . (0) + B.X.(O) + B.y.(0) 
1 1 J 2 ]' 3 4 M j 


= A j g 3 (z) + 2 A j B j g 3 (z)g 4 (z) + B j g 4 (z) 


(32) 


for j = 1...3, so that J is defined by 

3 N n2 f „ ... /A , A P z f 


1 VA *2 

; 'ho 


J U Lj k j 


0=1 


g 3 (z)g 4 (z) if 


iC 


*2 


2 . . dz 


(33) 


The integrals in (33) are 


*1 = 


'f -2Kz . 2. ,, 

e sm (az)dz 


w. 


'0 


P f -2Kz 
0 

= / sin(az)cos(az)dz 


( 34 ) 


0 ), 


0 


0 


J 


and 


(Eq. (34) continued on next page) 
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p f e -2Kz 2 
I~ = / cos axdz 

3 Jo 


Using the table of integrals, reference 135, these work out to be 


-2Kz -2Kz 

I i = ~Ti ft: — + — 5 TTT cos ( 2az ) " a sin ( 2az >) 

1 4Ka) 0 (4K I- 2 + 4aV 0 


_ -2kz / - K sin (2az) - a xos (2az) 

i 2 6 \ 22 

\ (4IC + 4a^)w 0 


I- = j s {cos az(-aK cos az + 2a sin az)} 

(4K Z + 4a ; Q 



